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ABSTRACT 

In  an  electrogasdynamic  (EGD)  generator  the  flow  of  a  neutral  medium 
carrying  unipolar  electrical  charges  constitutes  an  electrical  current. 
A  model  of  the  charge  flow  in  an  EGD  generator  was  constructed  for  use 
in  a  computer  program.  The  program  is  designed  to  solve  Poisson's  and 
Laplace's  equations  for  both  axi symmetric  and  two-dimensional  geometries. 
Schlieren  photographs  of  the  charge  cloud  were  used  to  determine  the 
charge  cloud  profile  required  by  the  program.  Computer  generated  pre- 
dictions agreed  with  three  known  solutions  to  Poisson's  equation. 
Computer  predictions  of  the  effects  of  space  charge  flow  modification 
were  obtained.  Space  charge  flow  was  modified  both  by  increasing  flow 
speed  and  by  manipulating  the  space  charge  electric  field.  Experimentally, 
this  was  accomplished  by  increasing  flow  stagnation  pressure  and  by 
application  of 'a  separate  controllable  electric  field.  Experimental 
results  compare  favorably  with  computer  predictions.  Some  measurements 
were  also  made  of  the  mobility  range  of  the  charged  particles. 
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I.  INTRODUCTION 

Electrostatic  phenomena  involve  many  everyday  experiences,  such  as 
walking  on  a  thick  rug  and  then  touching  a  grounded  object.  Practical 
applications  of  electrostatic  devices  range  from  Xerography  to  pollution 
control  devices  or  to  the  van  der  Graaf  generator  [1].  Recently,  interest 
in  electrical  power  generation  has  brought  attention  to  the  electrogas- 
dynamic  (EGD)  generator.  This  direct  energy  conversion  device  is  ana- 
logous to  the  van  der  Graaf  generator  and  operates  by  convecting  charges 
with  a  gas  against  an  electric  field  to  a  point  of  collection.  The 
charge  flow  constitutes  an  electrical  current.  The  charge  movement  is 
accomplished  by  an  exchange  of  momentum  between  the  flowing  dielectric 
medium  and  the  charges. 

The  common  denominator  of  electrostatic  devices  is  the  existence  of 
free  charges  in  a  neutral  medium  or  dielectric  --  for  the  present  pur- 
poses a  fluid.  The  flow  of  a  medium  with  unipolar  charges  is  termed 
space  charge  flow  and  the  understanding  of  this  flow  is  critical  to  the 
design  of  the  EGD  generator.  This  thesis  consists  of  a  computer  model 
of  the  space  charge  flow  to  be  used  as  an  EGD  generator  design  tool. 

The  computer  model  has  some  limitations  inherent  to  the  programming 
itself  and  these  are  discussed  in  Section  III.  In  addition,  parameters 
such  as  charge  mobility  and  velocity  as  well  as  the  inputs  of  geometry 
and  jet  configuration  are  needed.  These  are  obtained  with  an  experimental 
set-up  discussed  in  Section  IV.  The  program  is  used  to  predict  regions 
of  breakdown,  given  an  electrode  geometry  and  a  space  charge  cloud.  In 
addition,  the  program  is  used  to  predict  the  effect  of  guard  electrodes. 
These  electrodes  are  placed  for  the  purpose  of  'smoothing'  the  electric 
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field  in  the  conversion  region,  and  both  computer  predictions  and 
experimental  results  are  given  in  Section  VI  for  the  effect  of  the 
guard  ring. 


II.  THEORETICAL  CONSIDERATIONS 

The  EGD  process  requires  producing,  transporting  and  collecting  the 
charges  which  comprise  the  electrical  current.  One  means  by  which 
charges  may  be  produced  is  by  a  corona  discharge.  Basically  this  scheme 
consists  of  two  electrodes,  which  may  be  a  needle  and  a  concentric  ring; 
the  concentric  ring  or  attractor  electrode,  when  charged  to  a  high  posi- 
tive electrical  potential  with  respect  to  the  grounded  needle,  induces  the 
corona  discharge.  The  corona  breaks  down  a  neutral  gas  initially  into 
its  electrical  components:  molecular  ions  and  free  electrons.  The  corona 
needle  will  attract  the  positive  ions,  while  negative  charges  will  tend 
to  migrate  to  the  attractor  ring.  The  EGD  objective,  instead,  is  to 
transport  a  current  of  negative  charges  away  from  this  attractor  ring. 

An  ion's  susceptibility  to  movement  in  an  electric  field  is  termed 
it's  'mobility',  (k).  Thus  the  drift  velocity  of  the  ion  in  the  corona 
field  is  determined  by  the  product  of  the  ion  mobility  and  the  electric 
field  strength  to  which  the  ion  is  exposed.  Since  the  mobilities  of 
molecular  ions  and  free  electrons  are  high,  they  would  migrate  rapidly 
to  the  corona  electrodes  and  be  lost  to  EGD  power  generation.  However, 
when  a  supersaturated  vapor  is  present,  negatively  charged  droplet 
(micron)  sized  particles  are  nucleated  which,  for  a  given  charge,  will 
have  much  lower  mobility  than  an  ion  [2,5].  Hence,  the  charged  particle 
is  more  suited  to  EGD  applications  [6].  To  move  the  charged  particle 
against  an  electric  field,  the  particle  must  be  well  coupled  to  the 
dielectric  flow.  When  the  dielectric  gas  transports  the  charge  to  a 
downstream  collector,  against  an  electric  field  (E),  the  total  charge 
velocity  (v  )  may  be  described  as  the  sum  of  the  dielectric  velocity 
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(v  )  and  the  drift  velocity  (kE). 

Vp  =  v„  +  kE  (1) 

In  order  to  obtain  the  highest  degree  of  coupling,  the  drift  velocity 
of  the  charged  particle  must  be  low  compared  to  the  dielectric  velocity 
[4].  The  charged  particles  are  carried  out  of  the  charge  production 
section  against  an  electric  field  in  the  conversion  section,  to  the  charge 
collector.  This  collector  consists  of  an  equipotential  metallic  surface. 
A  path  through  a  load  to  ground  is  provided  to  complete  the  circuit 
started  at  the  grounded  corona  electrode.  See  Figure  (1)  for  a 
schematic  of  the  EGD  generator. 

There  are  several  limitations  that  pertain  to  the  conversion  process. 
For  an  EGD  generator  to  operate  as  intended,  the  neutral  gas  must  remain 
non-conducting.  This  implies  that  the  electric  fields  to  which  the 
dielectric  is  subjected  must  remain  below  the  dielectric  strength,  or 
breakdown  field  strength,  (E.  ),  of  the  gas.  Fields  above  the  dielectric 
strength  will  produce  regions  of  high  conductivity  through  ionization 
[3],  The  dielectric  gas  or  fluid  is  subjected  to  the  field  resulting 
from  the  applied  potential  through  which  the  charged  particles  are  carried 
and  to  internal  fields  generated  by  the  presence  of  the  charged  particles. 
Breakdown  from  external  fields  occurs  as  a  result  of  the  electrode  geo- 
metry. For  a  given  voltage,  the  existence  of  a  small  radius  of  curva- 
ture will  produce  a  large  electric  field.  The  field  will  exceed  the 
dielectric  strength  when  the  radius  of  curvature  is  small  enough.  The 
range  of  allowable  operating  voltages  is  thus  a  function  of  the  elctrode 
spacing.  In  turn,  the  spacing  of  the  electrodes  is  also  a  function  of 
the  fluid  dynamics  of  the  dielectric,  of  possible  charge  depletion 
mechanisms,  of  the  operating  pressures,  and  of  generator  size  requirements. 
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Total  charge  per  particle  is  limited  by  the  balance  between  the 
particle  surface  tension  and  charge  repulsion  forces,  and  is  known  as 
the  Rayleigh  limit.  In  practice,  however,  the  total  charge  per  particle 
is  governed  by  the  charge  production  method.  For  particles  charged  in 
a  corona  field,  and  which  grow  to  micron  size,  the  number  of  charges  is 
approximately  two  orders  of  magnitude  less  than  that  predicted  by  the 
Rayleigh  limit  [6].  These  limits  must  be  considered  in  the  modeling  of 
the  space  charge  flow  in  an  EGD  device. 

In  general,  the  equations  needed  to  describe  the  EGD  flow  will  include 
the  momentum,  continuity  and  energy  equations.  These  would  be  necessary 
to  determine  the  dielectric  jet  profile  and  hence  the  space  charge 
density  distribution.  However,  if  the  profile  can  be  determined 
experimentally,  and  if  the  energy  drained  from  the  flow  is  low,  the 
required  equations  are  reduced  to  three: 

Faraday's  Law 

r=  electric  field  (V/m)  « 

B  =  magnetic  flux  density  vector  (Wb/m  ) 


2 
p  =  space  charge  density  (C/m  ) 

e  =  dielectric  constant  for  vacuum  and  gas 

°   (V-m/A-s) 


V  •  Pg  Vp  =  0  (2) 

where  the  total  velocity  of  the  charged  particle  is  given  by: 

7  =  V  +  k(F^^„  +  E^^)  (3) 

p   <»     app   sc 

and  is  a  function  of  the  dielectric  flow  speed,  particle  mobility, 
applied  electric  field  (E^gnn^'  ^""^  ^^^  space-charge-produced  electric 
field  (E^^). 
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V  X  E  = 

.  3B 
at 

Gauss' 

Law 

V  .  E  = 

"e  ■ 

Charge 

Conservation 

For  the  EGD  generator,  B"  =  0,  and  Faraday's  Law  takes  the  form, 
V  X  "E  =  0.     This  allows  T  to  be  defined  as  the  negative  gradient  of  some 
potential   function     (}>,  so  that  E  =  -  V4).     Substituting  into  Gauss'   Law, 
the  resulting  equation  is  known  as  Poisson's  equation, 

A  -  -^  (4) 

This  equation  is  valid  in  the  two  regions  of  interest  within  the  EGD 
channel.  One,  the  neutral  gas  jet,  contains  the  charged  particles,  while 
the  second,  outside  the  gas  jet,  contains  no  space  charged.  In  the 
space-charge-free  region,  pp  =  0,  and  Eq.  (4)  reduces  to  the  more  famil- 
iar Laplace  equation.  Figure  (1)  depicts  the  two  regions. 

Although  the  kinetic  energy  of  the  injected  gas  jet  is  generally 
augmented  by  the  transfer  of  momentum  from  a  primary  flow,  this  method 
of  flow  augmentation  is  not  to  be  used  here  because  of  the  resulting 
unsteady  motion  of  the  jet  [7].  Rather,  the  injected  flow  is  augmented 
as  described  in  Section  IV. 

Reference  [8]  states  that  the  jet  profile  may  be  determined  by  con- 
sidering the  equation  of  conservation  of  charges.  The  radius  of  the  gas 
jet  (r)  may  be  described  as  a  function  of  the  initial  radius  (r  ), 
initial  charge  density  (Ppq) >  the  applied  electric  field  and  the  axial 


coordinate   (z). 

(5) 


-^o^f^F 


(v  +  kE^^  \ 
\  "    appj 


This  is  based  on  the  assumption  that  the  charge  density  is  a  function  only 
of  axial  distance  from  the  point  of  injection  and  that  the  particle 
velocity  due  to  coupling  with  the  neutral  gas  and  applied  field  is  much 
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greater  than  that  due  to  the  space  chrage  field,  i.e., 

V     +  k  E^,^»  k  F  ^  (6) 

»  app  sc 

Neglecting  space  charge  effects  disallows  any  radial  velocity.     Yet 

there  must  be  some  component  of  velocity.     For  flow  along  a  streamline, 

Crocco's  theorem  for  a  steady  state,  which  relates  enthalpy  (h   )  and 

entropy  (S)   gradients  to  vorticity,  is  given  by: 

v"    X  V  X  7    =  vh^  -  TvS 

.  .  00  CO  0 

where  v  x  7    =  w    is  the  vorticity.     For  the  streamline  defining  the 
boundary  of  the  gas  jet,  such  gradients  will  be  yery  large,  giving  rise 
to  vorticity.     The  radial   velocity  component  introduced  will   give  a 
larger  jet  profile  than  that  predicted  by  Eq.    (5),  leading  to  a  reduced 
charge  density  at  any  given  axial   distance. 

In  general ,  Equations   (2)  and  (4)  are  coupled  through  the  p     term, 
and  must  be  solved  simultaneously.     Expanding  the  charge  conservation 
equation: 

P  V 

T^If^P  ""T^lFPe  "^  Pelr^p  ""  ^p  If  ^e  =  ° 

For  particles  of  low  mobility,  and  fields  below  breakdown,  v  will 
predominate,  and  Eq.  (6)  will  hold.  If  the  jet  is  assumed  to  be  flowing 
normal  to  two  plane  parallel  electrodes: 

Then  an  assumption  of  an  initial  homogeneous  distribution  of  charge  at 
the  jet  entry  plane  will  remain  true  at  any  position  downstream.  Thus, 

3    -  n 
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and  the  charge  conservation  equation  becomes: 

V   '    io   w   )   =  0     —  V  +  V  T —  0   =0 

This  allows  the  cloud  shape  and  hence  the  charge  density  to  be  determined 
solely  by  the  fluid  dynamics  of  the  dielectric  jet,  and  uncouples  the 
charge  conservation  equation  from  Poisson's  equation.  Another  approach 
to  this  problem  is  to  determine  the  jet  profile  experimentally.  An 
approximate  analytical  description  is  obtained  by  fitting  a  parabola  to 
this  profile.  This  approach  is  convenient  in  this  work  and  the 
approximation  becomes  an  integral  part  of  the  computer  program. 

We  can  now  discuss  the  assumption  of  Eq.  (6)  that  the  charged  par- 
ticles are  affected  only  by  the  gas  jet  and  the  applied  field.  When 
charged  particles  are  injected  into  the  conversion  region,  they  will  face 
the  field  produced  by  the  charged  particles  already  in  the  conversion 
section.  Eventually,  these  charged  particles  will  face  fields  which  may 
repel  all  charges  but  those  with  the  highest  degree  of  viscous  coupling 
to  the  gas  flow.  In  this  manner,  space  charge  effects  provide  a  velocity 
"filter"  through  which  the  gas  jet  must  move  [9]. 

The  expression  for  current  density  (J)  in  the  dielectric  jet  is 
given  by: 

^e  p 
This  expression  would  require  a  greater  charge  density  as  space  charge 
effects  slow  down  the  charged  particles  in  order  to  maintain  a  given 
current  density.  But  as  mentioned  previously  on  Page  11,  breakdown  of  the 
dielectric  will  occur  if  the  space  charge  electric  field  strength  is 
greater  than  the  dielectric  breakdown  strength.  This  tradeoff  limits 
the  performance  of  the  EGD  generator. 
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III.     PROGRAMMING  CONSIDERATIONS 

The  computer  program  was  written  for  the  numerical  solution  of  Eq. 
(4),  since  closed  form  solutions  to  Poisson's  equation  are  difficult  to 
obtain,  except  in  very  simple  cases.     The  program  is  capable  of  solving 
Poisson's  equation  for  either  an  axisymmetric  or  two-dimensional   geometry. 

To  facilitate  the  solution,  Poisson's  equation  was  normalized  as 
described  in  Appendix  A.     In  the  axisymmetric  form,  the  equation  to  be 
solved  becomes: 

IAr  +  A^^  +  A^^^-Cp  (7) 

where  A  is  the  normalized  potential  and  R  is  the  normalized  radius.     The 
term  on  the  right  side  is  a  result  of  the  normalization  described  in 
Appendix  A.     In  the  two-dimensional   equation,  there  is  a  change  of 
variables  and  the  radius  term  drops  out  to  give  the  form: 

'^XX  "^  ^YY  "  "  ^^ 
By  covering  the  geometry  to  be  studied  with  a  uniform  square  mesh, 
the  problem  may  be  defined  in  terms  of  a  finite  number  of  discrete  points, 
at  each  of  which  the  solution  to  Eq.   (7)   is  desired.     The  boundaries  of 
the  mesh  are  determined  either  by  the  presence  of  electrodes  at  a  con- 
stant potential,  or  by  a  zero  normal   derivative  of  the  potential   cj). 
The  zero  normal   derivative  is  found  on  an  axis  of  symmetry,  on  a  boundary 
with  no  normal  current  flow,  such  as  an  insulator,  or  on  a  boundary 
between  two  electrodes  in  the  absence  of  space  charge  and  end  effects. 
End  effects  result  in  a  distortion  of  the  distributions  about  the  end  of 
a  pseudo-infinite  plate  because  of  the  radius  of  curvature  of  the  end. 
Figure   (2)  depicts  a  sample  geometry  with  appropriate  boundary  conditions. 
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The  derivation  of  the  finite  difference  equations  for  the  two  forms 
of  Poisson's  equation  is  given  in  Appendix  B.     After  solving  for  the 
potential   at  a  particular  point,  the  form  of  the  equation  is: 


^.j 


HL^Cp  +  A,^,    .  +A,    ,    .  +A.    .^,   +A.    .   T   +  ^A.    .^,   -A.    .J 
1+1, J       1-1, J       1  ,J+1       i,J-l       2R^   i,j+l       i,j-r 


where  the  subscripts  represent  mesh  row  and  column,  respectively.     For 
points  governed  by  a  boundary  condition  of  a  zero  normal   derivative,  as 
on  an  axis  of  symmetry,  the  general   equation  becomes: 

The  derivation  of  this  equation  is  also  given  in  Appendix  B. 

The  solution  method  chosen  is  an  iterative  process  known  as   "success- 
ive over-relaxation"   (SOR),  or  the  "extrapolated  Liebmann  method  [10]." 
Over-relaxation  introduces  a  factor  (w) ,  by  which  the  solution  at  the 
point  in  question,  as  a  result  of  the  previous  iteration,  is  taken  into 
account.     By  assigning  the  proper  value  to  this  factor,  convergence  rate 
of  the  solution  may  be  maximized;  however,  w  is   limited  by  the  following 

1   <  0)  <  2 

The  method  of  selecting  a  value  for  co  may  be  found  in  textbooks  covering 
iterative  solution  methods  such  as  Reference  [10].     For  this   thesis  w  is 
a  constant  equal   to  1.5.     The  form  of  the  resulting  equation  then 
becomes : 


a""^^  =  a"       +  ^ 


HL^Cp  +  a"  ,    .  +A'?"']    .  +  a"    .   T  +  a"-"^   . 

1+1, J       1-1, J         i,J-l  i,j+l 


^  2R\^,j+l       ^,j-l 


■4  a" 


1  ,J 
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This  shows  that  the  solution  at  a  point,  on  the  (n+l)st  iteration  is 

dependent  on  the  solutions  at  surrounding  points  obtained  on  either  the 

(n)th  or  (n+l)st  iteration.  Convergence  is  determined  by  comparing  a 

specified  value,  DELU,  to  the  maximum  potential  change  per  iteration. 

The  value  of  DELU  is  a  data  item  required  by  the  computer  program  and 

remains  constant  until  changed  by  the  user. 

This  simplification  of  the  governing  equations  introduced  in 

Section  II  allows  the  space  charge  density  to  be  described  as  a  function 

of  the  initial  space  charge,  of  jet  radius  and  of  the  initial  jet  radius, 

The  expression,  as  shown  by  Minardi,  is: 

r  \2 


(W 


^e   ^eo 

This  leads  to  a  charge  density  varying  with  the  longitudinal  coordinate 
only  and  a  constant  distribution  in  the  plane  normal  to  the  dielectric 
flow.  The  use  of  a  known  nozzle  radius,  i.e.,  intial  jet  radius,  of  a 
gas  flow  speed,  and  of  a  corona  current  (Ij),  allows  determination  of 
the  injected  space  charge  density. 

Also  to  be  considered  when  modeling  the  charge  distribution  is  the 
collection  process.  Charged  particles  cannot  continue  to  be  collected 
indefinitely.  Eventually,  at  some  point  downstream  from  the  initial 
point  of  collection,  all  but  a  negligible  part  of  the  charged  particles 
that  are  to  be  collected,  will  have  been  collected.  Further  charge 
contributions  to  the  collector  current  will  have  a  negligible  effect  on 
the  potential  distribution.  It  is. necessary  to  determine  where  this 
takes  place,  as  the  boundary  conditions  on  the  problem  include  a  zero 
charge  density.  Hence,  the  space  charge  contained  in  the  dielectric  jet 
is  not  to  be  carried  up  to  the  boundary  of  the  region  of  interest. 
Using  various  collection  lengths  in  the  computer  program,  it  was  found 
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that  allowing  space  charge  to  exist  over  80%  of  the  collector  length 
gave  results  which  changed  little,  less  than  3%,  as  this  collection 
length  was  increased. 

In  general,  the  jet  profile  does  not  coincide  with  a  mesh  point  at 
any  given  axial  coordinate.  However,  the  program  is  not  capable  of 
handling  a  variable  mesh  size.  In  order  to  differentiate  clearly 
between  the  regions  of  space  charge  and  no  space  charge,  the  program 
approximates  the  gas  jet  with  a  series  of  steps  in  order  to  meet  the 
constant  mesh  interval  limitation.  This  is  explained  with  the  help  of 
the  figure  below: 


r+Ar 


i,J 


HL 


i ,  J+1 


r+Ar 


i,J 


HI 
If  Ar  >  2"  >  the  gas  jet  is  considered  to  be  passing  through  the  point 

HI 
(i»j+l)>  and  if  Ar  <  2—,  through  the  point  (i,j).  The  charge  density  is 

computed  according  to  Eq.  (8),  giving  a  value  of  charge  density  which  is 
then  assigned  to  each  mesh  point  out  to,  and  including,  the  point  on  the 
approximated  jet  boundary.  This  introduces  some  error  into  the  calcula- 
tion process.  The  approximation  scheme  limits  the  displacement  of  the  ■ 
dielectric  jet  boundary  to  a  maximum  of  HL/2.  Then  the  maximum  error 

(°n.^^)   is  given  by: 
max    ^     •^ 

"rax  =  "L/2R 
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This  would  imply  that  the  maximum  error  would  occur  at  the  minimum  jet 
radius,  the  injector  nozzle.     By  selecting  the  appropriate  mesh  interval, 
the  radius  at  the  nozzle  can  be  made  to  coincide  with  a  mesh  point.     The 
minimum  radius  at  which  the  maximum  displacement  will  occur  is   then 
downstream  from  the  nozzle,     a        =  5%  is  typical  of  the  various  mesh 
sizes  used. 

The  validity  and  accuracy  of  the  program  was  determined  by  compari- 
son with  the  known  solution  for  several   problems.     A  geometry  consisting 
of  finite  parallel   plane  electrodes,  approximating  infinite  parallel 
plane  electrodes,  with  a  gap  spacing  of  2  cm  and  a  potential  difference 
of  15  kV  was     used  to  determine  if  the  program  would  predict  the  exist- 
ence of  electrical  breakdown  across  the  gap  as  specified  in  [11]. 
Because  of  the  method  of  normalization  used,  as  described  in  Appendix  A, 
the  result  of  entering  this  configuration  into  the  program  is  a  potential 
distribution  which  varies  linearly  from  zero  to  one,  and  a  field  strength 
of  one  throughout  the  gap,  indicating  that  breakdown  of  the  air  gap  is 
predicted.     To  attain  this  solution,  it  was  found  that  electrode  end 
effects  had  to  be  minimized,  or  a  yery  distorted  distribution  would 
result.     A  pseudo-infinite  plane  electrode  minimum  length  to  air  gap 
ratio  of  one  was  required  to  minimize  the  distortion.     Figure  (3)   depicts 
the  correct  potential  distribution  calculated  by  the  program  once, the 
parallel  plane  electrodes  were  properly  modeled. 

Results  of  the  program  were  also  compared  with  two  results  previously 
published  by  Minardi    [2,8],  concerning  space  charge  flow  between  parallel 
infinite  plane  electrodes;  values  of  the  input  data  used  were  the  same. 
A  comparison  of  the  graphical  output  appeared  to  be  equivalent  for  the 
two  geometries  compared. 
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IV.  DESCRIPTION  OF  APPARATUS 

The  need  to  input  certain  parameters  to  the  program  dictated  the  type 
and  extent  of  the  experimental  work  to  be  accomplished.  This  work  also 
had  to  be  relevant  to  the  EGD  facility  that  was  under  operation  at  NPS 
[4,5].  The  program  requires  both  the  electrode  geometry  of  interest  and 
the  dielectric  gas  jet  profile,  as  inputs,  so  that  charge  density  distri- 
bution may  be  calculated.  The  gas  jet  profile  is  also  the  means  of 
modeling  the  change  in  a  variable  such  as  the  dielectric  flow  speed.  The 
effect  of  an  increase  in  7  was  of  particular  interest  because  a  larger 

CO  ' 

V  would  increase  the  dominance  of  the  first  term  on  the  right  of  Eq.  (1). 

00 

Also  of  interest  was  the  effect  of  the  guard  ring  electrode  on  the  space 
charge  flow.  These  two  experiments,  increasing  7^  and  the  use  of  the 
guard  ring,  were  also  required  for  correlation  of  the  trends  in  experimental 
results  with  those  trends  predicted  by  the  computer  program. 

The  EGD  channel  to  be  modeled  consists  of  a  charged  aerosol  injector 
and  a  collector  probe.  See  Figure  (4)  for  dimensions.  The  vapor  for  the 
aerosol  is  generated  by  a  modified  Hotshot  model  MB-31  electric  steam 
boiler,  pressure  fed  with  distilled  water.  To  obtain  best  results,  it 
was  found  necessary  to  ensure  that  the  entire  steam  generating  system 
was  free  of  any  impurities.  The  corona  needle  is  grounded  through  a 
Simpson  microammeter,  while  the  attractor  ring  is  powered  by  a  Sorenson 
high  voltage  supply,  with  a  range  of  zero  to  ten  kV.  The  power  supply 
output  is  monitored  by  a  Sensitive  Research  high  impedance  voltmeter 
connected  across  the  corona  attractor  ring  gap.  The  collector  needle 
is  grounded  through  a  Simpson  microammeter  in  series  with  a  diode  to 
prevent  reverse  current  flow  [7]. 
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In  order  to  determine  the  jet  profile,  Schlieren  photographs  of  the 
steam  jet  were  taken.     The  setup  consists  of  mirrors,  a  collimated  light 
source,  a  knife-edge  and  a  camera.     Photographs  of  the  jet  were  taken 
with  stagnation  pressures  of  8,  10  and  14  psig.     A  horizontal   knife-edge 
was  used  to  more  clearly  define  the  vertical   gradient  across  the  jet 
boundary.     The  Schlieren  photographs  of  the  8  and  14  psig  jet  are  shown 
in  Figure  (5).     The  photographs  show  a  high  pressure  jet  (14  psig) 
approximately  100%  larger  than  the  given  by  Eq.   (5).     The  resulting  jet 
profiles  are  shown  in  Figure  (6). 

In  order  to  determine  the  effects  of  manipulating  the  electric  field 
distribution,  a  field  entirely  separate  from  that  of  the  corona  was 
applied  by  means  of  a   'guard  electrode',  in  the  form  of  a  ring  of  stain- 
less steel.     Internal   diameters  used  ranged  from  3/8  of  an  inch  to  3/16 
of  an  inch.     The  ring  was  mounted  on  a  plexiglass  insulator,  supported 
by  a  stainless  steel   rod.     The  rod  was  mounted  on  a  vernier  allowing  the 
position  of  the  ring  to  be  varied  along  the  longitudinal  axis  of  the 
collector  needle.     The  voltage  applied  to  the  guard  ring  was  supplied 
by  a  Spell  man  high  voltage  power  supply.     The  output  voltage  of  the  power 
supply  was  determined  by  calibration  of  the  power  supply  voltmeter  with 
the  Sensitive  Research  voltmeter.     This  was  necessary  because  the 
Spellman  voltmeter  does  not  indicate  line  voltage.     The  guard  ring 
schematic  is  shown  in  Figure  (7). 
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V.     EXPERIMENTAL  PROCEDURE 

Initial  experiments  with  the  guard  ring  were  carried  out  at  eight 
psig.     The  temperatures  used  corresponded  to  vapor  states  at  either 
side  of  saturated  vapor,  similar  to  previous   EGD  work  [4,5,7].     At  this 
low  pressure,  a  severe  deterioration  of  the   'convection  current'   ratio, 
(collector  current/corona  current),  was  encountered  as  ring  voltage  was 
increased.     This  appeared  to  be  the  result  of  too     low  a  dielectric 
velocity  and  hence  a  high  slip  parameter  (i.e.,  a  high  ratio  of  the  drift 
velocity  to  the  dielectric  velocity).     The  charged  particles  were  diverted 
from  the  collector,  being  attracted  to  the  guard  ring.     For  this  reason 
and  for  the  desire  to  improve  the  generator  performance  without  the  ring, 
it  was  decided  to  increase  the  flow  rate  of  the  dielectric  by  increasing 
the  reservoir  pressure.     Stagnation  conditions  were  changed  by  the 
introduction  of  a  separate   'dry'   gas  into  the  steam  boiler. 

The  baseline  characteristics  of  the  EGD  channel  were  determined  by 
supplying  steam  to  the  corona  discharge  at  eight  psig  and  240-245°F. 
Corona  voltage  was  varied  over  a  range  of  zero  to  three  kV,  and  values 
of  corona  current  and  collector  current  were  recorded.     The  collector 
probe  was  positioned  5  mm  downstream  from  the  nozzle  exit  plane.     Dry 
gases  used  were  both  Np  and  A.     The  choice  of  gas  appeared  to  make  no 
difference  in  the  results.     Introducing  the  dry  gases  was  done  by 
regulating  the  flow  of  gas  from  a  standard  high  pressure  gas  bottle. 
The  gas  was  introduced  into  the  boiler  through  acheck  valve  operating 
at  10  psig.     With  the  boiler  operating  at  8  psig,  the  high  pressure  gas 
was  used  to  raise  the  operating  pressure  of  the  boiler  to  approximately 
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13.5  psig.  After  allowing  the  rate  of  steam  production  to  stabilize, 
resultant  stagnation  pressure  was  14  psig. 

To  determine  the  effects  of  the  'guard  electrode',  a  baseline 
performance  with  no  guard  ring  was  run  at  14  psig.  Once  this  was 
established,  the  trailing  edge  of  the  guard  ring  was  placed  in  one  of 
three  positions,  either  2  mm  upstream,  coincident  with,  or  2  mm  down- 
stream of  the  collector  needle  tip.  At  each  guard  ring  position,  the 
applied  potential  was  varied  over  a  range  of  0  -  3.5  kV,  and  effects  on 
collector  current  were  recorded.  Rings  of  various  sizes  were  used. 
When  rings  smaller  than  3/8  of  an  inch  diameter  were  used,  breakdown 
occurred,  apparently  caused  by  stray  droplets  providing  a  direct  path 
for  current  transmission.  Condensation  of  the  steam  inside  the  ring 
was  also  observed,  and  was  a  possible  factor  contributing  to  breakdown. 
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VI.     RESULTS  AND  DISCUSSION 

The  introduction  of  a  high  pressure  gas  into  the  system  resulted  in 
collected  currents  of  two  to  three  times  greater  than  the  low  pressure 
steam.     Figure  (8)   shows  the  baseline  (8  psig)  corona  current-voltage 
curve  and  the  corona  current-voltage  curve  for  the  high  pressure  steam 
jet.     The  two  curves  are  representative  of  the  range  of  corona  operation 
Figure  (8)  also  shows  the  variation  of  the  value  of  the   'convection 
current'   ratio. 

In  the  expression: 

V „  =  V    +  k  r  ^^  +  k  El^ 
p        «>  app  sc 

the  term  7^  may  be  given  by  the  following  expression  from  Shapiro  [12]. 
where 


p 

,.-! 

V? 

M  =  - 

v" 

00 

For  the  pressures  considered,  8  psig  and  14  psig,  the  respective  Mach 
numbers  are  0.832  and  1.02.     These  correspond  to  flow  speeds,  at  the 
injector  nozzle  exit  plane,  of  402  m/s  at  P     =  22.7  psia,  and  483  m/s 
at  P     =  28.7  psia.     From  the  Schlieren  photographs  the  profile  of  the 
jet  was  determined.     Assuming  continuity,  this   leads   to  a  prediction  of 
the  variation  of  flow  speed  in  the  flow  direction. 

The  increased  effectiveness  of  the  high  pressure  jet  may  be  shown 
by  recalling  the  slip  parameter,  6,  defined  to  be  the  ratio  of  the  par- 
ticle drift  velocity,  kF,  to  the  steam  jet  flow  speed.     In  terms  of  the 
non-dimensional   electric  field  output  of  the  computer  program,  E, 
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the  slip  parameter  is 


a=!^ 


For  a  jet  with  particles  well  coupled  to  the  flow,  6  will  be  less  than 
one.  To  determine  6,  the  value  of  E  may  be  solved  for  by  the  computer 
program.     From  Eq.   (8),  the  value  of  p^^,  the  initial   charge  density,  may 

be  determined.     For  the  low  and  high  pressure  jets,  this  results   in 

3 
0.0102  and  0.0084  C/m  ,  respectively.     Using  this  density  as  the  program 

input  parameter  RHOZRO,  an    |E|   field  distribution  may  be  obtained.     This 
was  done  for  both  low  and  high  pressure  cases.     Other  significant  para- 
meters used  were  attractor  voltage  of  2300  volts,  and  breakdown  field 
strength  of  3  x  10     V/m.     Figures   (9)  and  (10)   are  plots  of  the  equi- 
potential  surfaces  resulting  from  the  computer  solution  of  the  low  and 
high  pressure  jets.     Using  these  results,  the  variation  of  6  along  the 

centerline  of  the  EGD  channel  was  obtained,  and  is  plotted  in  Figure  (11) 

-5     2 
The  values  of  6  are  based  on  an  assumed  mobility  of  4  x  10       m  /V-s, 

a  value     representative  of  the  range  of  mobilities  of  the  EGD  charged 

particle  [7]. 

Figure  (11)  shows  that  everywhere  along  the  centerline  of  the  EGD 

channel  ,  with  the  exception  of  a  single  point,  the  high  pressure  jet  has 

a  lower  slip  than  the  low  pressure  jet.  The  single  point  that  differs 

may  well  be  the  result  of  the  inaccuracies  of  the  numerical  solution. 

The  use  of  a  higher  value  of  mobility  would  tend  to  give  values  of  6 

greater  than  one.  However,  if  the  conjectured  distribution  of  [5]  is 

approximately  correct,  particles  of  higher  mobility  are  a  small  part  of 

the  total  number  of  charged  particles. 
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From  Eq.  (3),  the  total  E  field  has  two  contributions.  One,  E   , 
is  due  entirely  to  the  geometry  and  voltages  applied  to  various  parts 
of  the  EGD  geometry.  This  can  be  changed  only  by  a  physical  change  in 
the  geometry  or  applied  voltages.  The  potential  and  F  field  distribution 
resulting  from  the  geometry  alone  may  be  solved  for  by  Laplace's  equation, 
The  Laplacian  solution  represents  the  potential  or  E  field  distribution 
unaffected  by  the  presence  of  space  charge.  Figure  (10)  is  the  solution 
to  Poisson's  equation,  and  it  shows  the  added  distortion  of  the  dis- 
tribution due  to  the  presence  of  space  charge,  when  compared  to  the 
Laplacian  solution.  Figure  (12).  The  added  distortion  appears  in  Eq. 
(3)  as  the  second  contribution,  F  ,  the  space  charge  field.  Figure  (13) 
depicts  the  |E|  field  distribution  corresponding  to  the  potential 
distribution  of  Figure  (10).  Since  the  space  charge  is  of  the  same 
polarity  as  the  charged  particles,  their  movement  is  opposed.  Elimina- 
tion or  reduction  of  F  would  enhance  charged  particle  velocity, 
making  the  collection  easier.  Application  of  a  potential  to  the  guard 
ring  is  designed  to  impose  an  additional  external  field  on  the  flow. 
However,  if  the  field  overcompensates  for  the  space  charge  field,  and 
the  ring  becomes  a  dominant  sink  for  the  charged  particles,  the  purpose 
of  the  ring  is  defeated. 

Figure  (14)  shows  the  increase  in  collector  current  as  a  result  of 
using  the  guard  ring.  The  results  are  shown  over  a  range  of  voltages 
at  the  three  positions  used.  It  appears  that  the  position  achieving 
the  greatest  increase  is  that  of  the  trailing  edge  plane  coincident  with 
the  tip  of  the  collector  needle.  From  Figure  (10),  this  is  also  the 
area  of  greatest  space  charge  distortion,  implying  that  the  effects  of 
the  ring  are  concentrated  in  a  small  area. 
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Figures  (15),  (16)  and  (17)  are  graphical  computer  outputs  of 
potential  distribution  for  a  grounded  ring,  and  potentials  of  1  and  2  kV. 
The  experimental  results  also  show  the  proper  trend  in  collector  current 
as  ring  voltage  is  increased.  That  is,  the  current  slowly  increases  to 
a  broad  maximum,  and  then  drops  off  with  further  voltage  increases,  the 
maximum  occuring  at  about  1  kV. 

The  computer  predicted  results  of  using  the  guard  ring  is  shown  in 
Figure  (18)  as  a  plot  of  centerline  Laplacian  E  minus  the  Poisson  E 
field.  As  this  difference  approaches  zero,  the  particular  guard  ring 
distribution  approaches  the  optimum  Laplacian  distribution.  The  computer 
results  indicate  that  the  E  field  is  definitely  modified  by  the  presence 
of  the  guard  ring.  The  predicted  optimum  ring  voltage,  corresponding  to 
maximum  collector  current,  appears  to  lie  between  ground  and  1  kV.  The 
point  common  to  all  curves  Ties  at  the  midpoint  of  the  ring,  and  is  prob- 
ably due  to  the  uniformity  of  the  guard  ring  electric  field  at  this 
point. 
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VII.  CONCLUSIONS 

The  computer  program  appears  to  present  an  accurate  picture  of  the 
potential  and  electric  field  distributions  within  an  EGD  device.  This 
is  borne  out  by  the  computer  prediction  of  experimentally  observed  trends 
in  EGD  channel  performance.  The  'convection  current'  ratio,  Iq/I-i-j  was 
found  to  be  affected  by  the  steam  reservoir  conditions.  In  addition, 
computer  and  experimental  results  agree  and  show  that  increased  EGD 
performance  can  be  attained  by  suitable  manipulation  of  the  electric 
field  distribution,  and  hence  of  the  space  charge  flow. 

The  agreement  found  between  analytical  and  experimental  results 
suggests  that  the  use  of  the  computer  program  as  a  design  tool,  par- 
ticularly when  modeling  the  space  charge  flow  in  a  high  pressure  jet, 
is  feasible. 

To  be  completely  general  in  approach,  some  further  account  should 
be  taken  of  particle  mobility.  Since  the  drift  velocity  appears  to  be 
negligible  at  high  pressures,  this  is  especially  true  when  working  with 
particles  of  high  mobility,  and/or  dielectric  jets  at  low  speeds.  In 
this  case,  the  slip  parameter  could  go  well  above  1,  and  the  present 
assumption  of  zero  slip  will  give  completely  invalid  answers.  Some 
preliminary  work  on  determining  the  mobility  range  of  the  charged 
particles  used  in  these  experiments  is  shown  in  Appendix  D. 
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VIII.  RECOMMENDATIONS 

In  addition  to  the  measurement  of  mobility,  the  difference  equations 
should  be  modified  to  allow  for  a  variable  mesh  size.  Equations  of  this 
type  may  be  found  in  a  textbook  such  as  Reference  [10].  This  will  allow 
a  greater  generality  in  selecting  geometries  to  be  studied,  and  will  also 
do  away  with  a  great  deal  of  the  error  involved  in  discrete  jet  profile 
approximation.  It  will  almost  certainly  mean  greater  computer  storage 
requirements. 

As  a  design  tool,  the  program  can  be  used  to  determine  the  optimum 
guard  ring  geometry  and  voltage.  From  the  curves  of  Figure  (18),  the 
solutions  for  the  various  voltages  diverge  from  the  optimum  in  the 
vicinity  of  the  collector  needle.  A  guard  ring  geometry  capable  of 
affecting  this  divergence,  yet  maintaining  the  distribution  of  the 
remainder  of  the  1  kV  curve,  would  appear  to  be  an  improvement. 

It  should  also  be  possible  to  use  the  program  as  a  design  tool  to 
optimize  both  injector  nozzle  and  collector  geometries.  A  nozzle 
geometry  allowing  for  the  use  of  smaller  guard  ring  sizes  may  allow 
for  more  effective  guard  ring  performance,  yet  forestall  breakdown. 
An  optimum  collector  geometry  may  exist  which  will  allow  for  greater 
charge  collection  efficiency  resulting  simply  from  the  geometry,  rather 
than  by  flow  manipulation  to  achieve  the  higher  collector  currents. 

Other  applications  outside  the  EGD  field  may  also  exist.  For 

instance,  the  program  may  be  used  to  study  the  region  of  space  charge 

known  as  an  electrode  'sheath'  which  divides  a  neutral  plasma  and  an 

electrode.  Also,  Reference  [4]  covers  a  study  of  EGD  control  of  separated 

flow.  This  study  may  be  extended  through  use  of  the  program  to  determine 

optimum  electrode  geometries  for  control . 

30 


APPENDIX  A:  NORMALIZATION  PROCEDURE 


The  vector  form  for  Poisson's  equation  is: 

V  d)  = 

After  expanding  in  cylindrical   coordinates: 

ill  +  O.   +  l_i   =     ^ 


By  defining: 


where: 

h     =     characteristic  channel   length  of  EGD  geometry 

V     =     maximum  potential 
and  substituting  into  the  various   terms: 

h(^)     h  9(-)  h  R  ai^        h^  8(^)      ^ 

n 

.2 
similarly  for  — %  .     By  combining  and  clearing  the  left  side: 

1  ^  % 

±A      +A        4.n         =_   o 

R  ^R       ^RR       ^ZZ  Ve 

0 

Defining:  « 

r     _  eo     .  _    ^e 

^     ~  ■  ~v7         '  P  ~ 


Ve         '  ^         p 

0  ^60 

The  normalized  equation  becomes: 
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Similarly  for  the  two-dimensional   form  of  Poisson's  equation.     By 
defining: 

X  =  —  •  Y  =  ^ 

The  two  dimensional  form  becomes:  "*" 

The  electric  field  may  be  normalized  by  defining: 

where:     E  =  the  normalized  value  of  E 

E.  =     the  dielectric  breakdown  field  strength 

from:  E  =  -v^ 

^^  -   S^  =  -  vA 


V      V 
then:  E  =  -  r^  vA 
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APPENDIX  B:     FINITE  DIFFERENCE  EQUATION  DERIVATION 


J 


■^  I 


i-l 
J+1 


i-l 
J 


i-l 
J-1 


i  J+1 


i,J 


HL 


iTT 


i+1 
j+1 


i+1 
J 


i+1 
J-1 


Knowing  the  value  at  any  point  in  the  mesh,  the  value  at  another 
point,  nearby,  may  be  found  by  means  of  a  Taylor  series.     By  writing  the 
series  in  the  J  direction,  expressions  for  A.    .-,   and  A.    .   -,  may  be 
found.     These  are: 


A  =  A        +  HI  A '       +  ^^     A ' '     +  A ' ' '   + 


A.  .  T  =  A.  .  -  hla:  ,  +  ^  a:  ' .  -  ^  a:  ".  +  . . . 


(1) 


(2) 


Solving  for  AI    .,  neglecting  higher  order  terms  where  primes  indicate 
the  derivative  with  respect  to  the  direction  of  interest: 


A* 


_  ^ ,  j+1 '  ^ ,  j-1 


To  determine  the  finite  difference  equation  for  the  second  partial 
derivative,  add  (1)  and  (2).     This  results   in: 

A,-    ^4-1  +  A.    .  1   =  2  A.    .  +  HL^  A!'  .  +  ... 
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Solving  for  A!'.,  again  neglecting  higher  order  terms 
1  »J 


a: 


i,j 


_  ^,J-H  ^  ^•  J-1  '  ^  ^iJ 


(4) 


By  writing  the  series  in  the  I  direction,  solving  for  values  of  A.^,    . 
and  A.   -.    .,  and  adding  the  resulting  equations,  the  second  derivative 
in  the  I  direction  may  be  determined.     The  resulting  equation  is: 


A' ' 


(5) 


By  combining  (3),   (4),  and  (5),  and  solving  for  A.    .,  a  finite  difference 
form  of  Poisson's  equation  in  axisymmetric  form  is  obtained.     The 
axi symmetric  form  is: 


^ .  j "  4 


HL' 


■^p"^.l,j"^-l,j^^,3.l^^,M 


^  2R  ^^.j+l "  ^,j-l^ 


(6) 


For  points  along  a  boundary  where  the  normal  derivative  is  equal  to 
zero,  the  equations  for  the  value  on  the  boundary  may  be  derived  in  the 
same  manner  and  cast  in  the  form  of  second  order  forward  or  backward 
difference  equations.  Assuming  (i,j-l)  is  on  the  boundary  of  interest, 
the  series  equations  become: 

A,-  n^i  =  A.  •  1  +  2  hla:  ,  T  +  ^^^  l\'.'  .  ,  +  ^.^  ^\'\   T  +  ...(7) 
1  ,J+1    1  ,J-1       i,J-l    2!    i,j-l    3!    i,j-l     ^  ' 


A.  .  =  A.  .  T  +  hla:  .  1  +  ^  a:  ' .  ,  +  ^  a:  ".  ,  +  ... 

1  ,J    i.J-1     i,J-l   2!   i,j-l   31   1  ,j-l 
Multiplying  (8)  by  4  and  subtracting  from  (7): 


(8) 
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Neglecting  higher  order  terms: 


N.J-1 


-  2  HL 


(9) 


Setting  the  normal   derivative  equal   to  zero,  and  solving  for  A.    •   -, : 


^.j-l   =  1-333  A,,  ^j-    .333  A.  ^j,, 


(10) 


The  successive  over-relaxation  method  (SOR)   uses  previously  calculated 
values  to  speed  up  convergence.     For  a  point  inside  the  mesh,  the  final 
difference  equation  becomes: 


a""^^  =  a"      +  ^ 


+  a"  ...  +  a"-"! 


Hl'cp-aJ,1,-A^!1,j    ■^i,j,i    -^i.j.i 


n+1 


+  nk    fA"  -  a"    '     )   -  4A" 

^  2R    ^^.j+1      ^,j-l^      ^^•,j 


(11) 


For  a  point  on  a  boundary  of  zero  normal  derivative,  the  equation 
becomes : 


fln+l   _  a"    4. 

^,j-l  ■^J-l  "■" 


1.333  A^^.-  .333  A^^.,^  -  A^^. 


(12) 
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APPENDIX  C:  PROGRAM  DESCRIPTION 

The  computer  program  is  designed  to  solve  a  finite  difference  form 
of  Poisson's  equation,  either  in  axisymmetric  or  two-dimensional  form, 
at  each  node  defined  by  a  uniform  square  mesh  covering  the  geometry  to 
be  studied.  In  order  to  do  this,  the  program  must  be  able  to  do  the 
following: 

1.  Determine  which  finite  difference  equation  is  to  be  solved. 
This  is  a  function  of  the  geometry  and  the  node  position  within 
the  mesh. 

2.  Determine  and  store  charge  distribution  at  each  node. 

3.  Solve  the  difference  equation  and  store  the  answer. 

4.  Determine  when  a  solution  to  the  problem  has  been  attained. 

5.  Provide  the  results  in  an  easily  interpreted  form. 

The  program  is  made  up  of  two  major  parts:  preparation  and  calcula- 
tion, and  output  of  results.  In  general  at  each  point  in  the  EGD  channel 
defined  by  the  uniform  mesh,  there  exists  a  charge  density  and  a  potential 
The  charge  density  may  be  zero.  This  can  be  the  case  inside  of  the  jet, 
as  well  as  outside.  The  voltage  may  either  be  the  result  of  an  electrode 
presence  or  a  potential  surface  between  electrodes.  If  the  voltage  is 
due  to  an  electrode  presence,  this  value  will  remain  constant,  otherwise 
it  will  vary  as  a  function  of  the  charge  density.  The  program  must  be 
able  to  make  this  differentiation.  In  all,  the  program  must  know  the 
potential  at  a  point,  the  charge  density,  and  whether  the  point  defines 
a  portion  of  an  electrode.  In  order  to  simplify  the  program  coding,  the 
program  stores  all  values  of  potential  in  one  array,  "A".  The  values  of 
charge  density  are  stored  in  a  separate  array,  "G",  of  the  same  size. 
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This  array  is  also  used  to  store  all  information  required  to  determine 
whether  the  potential  is  that  of  an  electrode,  the  potential  description 
value  (PDV),  or  is  to  be  changed  as  a  function  of  the  charge  density. 
This  can  be  done  since  a  point  of  constant  potential  will  be  an  electrode, 
and  a  charge  density  cannot  exist  at  the  same  point. 

Th.e  program  determines  convergence  of  the  solution  method  by  compar- 
ing answers  from  the  (n+l)st  iteration  with  those  obtained  on  the  (n)th 
iteration.  This  will  reveal  the  maximum  change  in  potential  as  a  result 
of  the  iteration.  If  the  maximum  is  less  than  some  designated  value,  the 
iterative  procedure  is  stopped.  This  method  requires  storage  for  the 
answers  generated  by  the  (n)th  iteration.  These  are  stored  in  array  "B". 
The  value  of  the  potential  resulting  from  the  previous  solution  is  required 
for  the  SOR  method.  This  value  is  supplied  from  array  "B".  Hence  the 
variables  required  for  the  solution  of  the  finite  difference  equation 
are  taken  from  three  separate  sources.  A  superposition  of  arrays  "A" 
and  "B"  describes  the  complete  situation  at  any  node  within  the  EGD 
channel . 

The  program  has  been  written  for  use  on  an  IBM  360/67  Computer 
System.  The  output  is  in  both  printed  and  plotted  form.  The  plotting 
was  accomplished  on  the  CALCOMP  Model  765.  The  program  consists  of 
the  following  subroutines: 

PURPOSE 


SUBROUTINE  NAME 
MAIN  (Main  Program) 

GSPRD 


GENFIL 
CALCl 


General  coordination  of  preparation 

and  calculation  phase  and  output 

phase. 

Calculates  cloud  radius.  Determines 

charge  density  at  each  point,  and 

stores  in  "G"  array. 

Reads  geometry  input. 

Calculates  potential.  Determines  E 

field  components  and  modulus. 
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WRITER  Prints  values  of  potential,  E  field 

and  E  field  components,  determined 
in  CALCl. 

CLVL  Determines  lines  of  constant  potential 

and  constant  E  field  to  be  plotted. 

FLOP  Orients  arrays  to  coincide  with  output 

labeling  convention. 

OUTPUT  General  coordination  of  output  phase. 

CONTUR  Coordination  of  labeling  and  plotting 

of  graph,  and  preparation  for  plotting 
results.  Modified  NPGS  Computer 
Library  subroutine,  which  allows  the 
superposition  of  three  sets  of  contours 
on  one  plotted  output. 

All  data,  except  the  electrode  geometry  is  read  in  under  NAMELIST 
control.  This  requires  a  special  format  for  all  cards,  but  provides  for 
a  more  flexible  means  of  handling  the  data  in  that  there  are  no  fixed 
fields  in  which  the  variables  must  appear,  nor  must  the  arrangement 
coincide  with  that  of  the  NAMELIST  statement. 

The  first  data  card  required  by  MAIN  is  the  card  denoting  the  appear- 
ance of  the  data  group.  It  has  the  format,  starting  in  column  2;  &PARAM. 
PARAM  is  the  NAMELIST  name.  The  first  column  of  all  data  cards  is  ignored. 
Each  data  item  may  be  separated  from  the  following  by  a  comma.  However, 
no  comma  must  appear  between  the  NAMELIST  name  and  first  data  item. 
More  than  one  data  item  may  appear  on  a  card.  Each  data  item  is  in  one 
of  the  two  forms  below. 

1.  "Variable  name"  =  FORTRAN  constant 

2.  "Array  name"  =  a  set  of  FORTRAN  constants,  separated  by  constants 
[13]. 

The  variables  read  in  are  listed  below.  The  first  five  are  explained 
with  the  help  of  the  representative  geometry  in  Figure  (  2). 
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NAME 
NATRAD 

NBLOCK 


NDLPNT 
NDLRAD 

NOZPOS 

BRKDWN 
CONVER 
DELU 


GAPPOT 


TYPE 
1*4 

1*4 


1*4 
1*4 

1*4 

R*4 
R*4 
R*4 


R*4 


IH 

1*4 

IW 

1*4 

NCOLS 
NDIMEN  =  0 

1*4 
1*4 

=  1 

NGLOSS  =  0 

1*4 
1*4 

DESCRIPTION 

Mesh  column  number  of  attractor  inner 
radius.  This  is  the  point  defined  as 
(6,3).  Then  the  value  assigned  to 
NATRAD  is  3. 

The  number  of  discrete  points,  lines, 
or  blocks  of  points  required  to 
completely  describe  the  charge  dis- 
tribution or  geometry  of  the  problem. 
For  the  figure  shown,  the  attractor 
ring  may  be  defined  by  a  block  of 
points,  a  line  and  a  single  point. 
The  corona  needle,  and  collector 
needle  may  each  be  defined  by  two 
lines  of  points.  In  addition,  if  it 
is  assumed  that  the  collector  needle 
collects  all  the  charge,  a  block  of 
points  can  be  used  to  define  the  area 
within  the  dielectric  jet  where  in  no 
charge  exists.  Thus,  for  the  geometry 
and  charge  distribution  of  the  problem 
shown,  NBLOCK  =  8. 
Mesh  row  position  of  corona  needle 
tip.  For  the  geometry  shown,  NDLPNT  = 
Mesh  column  number  of  corona  needle 
centerline.  For  geometry  shown, 
NDLRAD  =  1 . 

Mesh  row  number  of  attractor  exit 
plane.  For  the  geometry  shown, 
NOZPOS  =  6. 

Breakdown  field  strength  of  dielectric 
gas.  (volts/meter) 

Characteristic  channel  length  of  EGD 
geometry,  (centimeters) 
Value  compared  with  maximum  difference 
between  iterations.  If  the  maximum 
difference  is  less  than 
of  DELU,  convergence  of 
scheme  is  signaled. 

Maximum  fixed  potential  covered  by  the 
mesh.  In  general,  the  attractor  will 
be  at  the  maximum  potential. 
Desired  height  of  output  plot  measured 
along  longitudinal  coordinate  (IH  <  15) 
Desired  width  of  output  plot  measured 
along  radial  coordinate  (IW  <  9) 
Number  of  mesh  columns  (NCOLS  ^  61). 
Indicates  problem  is  a  two-dimensional 
geometry  in  constant  area  channel. 
Indicates  problem  is  axi symmetric. 
No  charge  density  decrease  accounted 
for  in  the  conversion  section  other 
than  by  charge  cloud  spreading. 


desired  value 
the  iterative 
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NGLOSS  =  1 

1*4 

NL 

1*4 

NLI 

1*4 

NOLINE  =  -1 

1*4 

=  0 

1*4 

=  1 

1*4 

NROWS 
RHOZRO 

1*4 
R*4 

Some  loss  mechanism  (recombination, 

etc.)  and  cloud  spreading  will   be 

accounted  for  in  determining  charge 

cloud  densities. 

Number  of  potential   contours  desired 

in  plotted  output  (NL  ^  20). 

Number  of  E  field  contours  desired  in 

plotted  output  (NLI^<  20). 

Only  potential   and  E  field  contours 

will   be  plotted. 

Only  physical   geometry  (fixed  potential) 

points  will   be  plotted. 

Physical   geometry  will  be  superimposed 

on  contour  plot. 

Number  of  mesh  rows   (NROWS  <  81) 

Charge  density  at  exit  plane  of 

attractor  (coulombs/cubic  meter) 

The  last  card  in  the  NAMELIST  signals  the  end  of  the  list.     The  format 
is  &END.     Again  the  symbol  ampersand  must  appear  in  column  2  of  the  data 
card. 

The  jet  radius  and  charge  density  distribution  are  calculated  by 
subroutine  GSPRD  according  to  the  charge  profile  given  by  statement 
numbered  35.     In  the  case  of  a  two-dimensional   geometry,  only  the  charge 
distribution  is  calculated,  as  the  jet  is  assumed  to  fill   the  entire 
two-dimensional   channel.     The  charge  density  will   remain  constant  through- 
out the  flow  in  the  channel,  unless  some  charge  loss  mechanism  is 
introduced,  signaled  by  the  value  of  NGLOSS.     The  jet  profile  and  charge 
distribution  for  a  geometry  with  an  off  center  corona  needle  is  also 
determined  by  GSPRD. 

Subroutine  GENFIL  is  used  to  input  the  specific  values  of  potential 
and  charge  density  or  PDV.     This  is  accomplished  with  the  help  of  the 
variable  NBLOCK,  which  is  used  as  a  DO-LOOP  parameter.     The  loop  causes 
the  program  to  read  the  exact  number  of  data  cards  defined  by  NBLOCK. 
Since  each  block  defined  by  NBLOCK  contains  points  of  common  value,  the 
entire  block  can  be  represented  by  a  single  data  card.     The  format  for 
these  data  cards  is  4I3,2F12.4 
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The  variables  read  in  by  GENFIL  are: 

NAME  TYPE      DESCRIPTION 

Nil  1*4      Mesh  row  nunter  of  beginning  of  line 

or  block. 
NI2  1*4      Mesh  row  number  of  end  of  line  or 

block. 
NJl  1*4      Mesh  column  number  of  beginning  of 

line  or  block. 
NJ2  1*4      Mesh  column  number  of  end  of  line 

or  block. 
AV  R*4      Normalized  potential  common  to  all 

points  encompassed  by  Nil,  NI2, 

NJl  and  NJ2. 
GV  R*4       EITHER:  Normalized  charge  density 

common  to  all  mesh  points  encompassed 

by  Nil,  NI2,  NJl  and  NJ2. 

OR:   Potential   description  value  (PDV) 

if  the  defined  points  are  at  a  fixed 

potential. 

An  arbitrary  value  of  99.  is  the  PDV  used  to  describe  a  point  of  fixed 
potential.  No  arithmetic  operations  will  be  carried  out  at  this  point. 
Note  also  that  Nil  <  NI2,  and  NJl  ^  NJ2,  since  these  variables  are  the 
beginning  and  end  of  a  DO-LOOP.  In  the  case  of  a  point.  Nil  =  NI2  and 
NJ  1  =  NJ2.  In  the  case  of  a  line  of  points,  either  Nil  =  NI2  or  NJl  = 
NJ2,  depending  on  the  orientation  of  the  line.  The  values  assigned  to 
the  above  variables  to  describe  the  representative  geometry  are  shown 
bel ow : 

Nil     NI2    NJl     NJ2    AV     GC    Corresponding  to 


1 

6 

5 

15 

1.0 

99. 

NBLOCK  1 

5 

6 

4 

4 

1.0 

99. 

NBLOCK  2 

6 

6 

3 

3 

1.0 

99. 

NBLOCK  3 

1 

4 

1 

1 

0.0 

99. 

NBLOCK  4 

1 

3 

2 

2 

0.0 

99. 

NBLOCK  5 

14 

20 

1 

1 

-0.1 

99. 

NBLOCK  6 

15 

20 

2 

2 

-0.1 

99. 

NBLOCK  7 

17 

20 

3 

6 

0.5 

0.0 

NBLOCK  8 

At  the  same  time  that  the  values  of  AV  and  GV  are  being  loaded  into  the 
appropriate  arrays,  the  coordinates  of  those  points  with  a  PDV  of  99. 
are  being  converted  into  plotter  pen  position  coordinates  to  aid  in 
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plotting  the  output  requested  by  NOLINE.  The  pen  position  coordinates 
are  stored  in  the  appropriate  vectors,  XLIN  and  YLIN  for  later  use  during 
the  output  phase. 

Subroutine  CALCl  solves  the  appropriate  form  of  the  finite  difference 
equations,  depending  on  the  geometry  variable  NDIMEN,  and  location  of  the 
point  in  question.  In  addition,  once  the  potential  distribution  has  been 
calculated,  CALCl  determines  the  E  field,  and  the  longitudinal  and  radial 
components  of  the  E  field.  The  value  of  the  E  field  and  its  components 
are  normalized  with  respect  to  the  breakdown  field  strength  of  the 
dielectric,  BRKDWN,  while  the  values  of  potential  are  normalized  with 
respect  to  the  maximum  fixed  potential,  GAPPOT.  The  derivation  of  the 
normalization  of  E  is  shown  in  Appendix  A,  which  gives: 

.   ^  =  -  ^  VA 

The  characteristic  EGD  geometry  length  'h'  is  read  into  the  program  by 
the  value  of  the  variable  CONVER.  Defining: 

$ALPH  =  '^ 


^h 


|E|     =  $ALPH    |-vA| 

The  value  of  the  constant  $ALPH  is  calculated  by  MAIN,  and  is  passed  to 
CALCl   as  a  calling  argument.     By  normalizing  E  in  this  manner,  a  quick 
glance  at  the  results  will   show  whether  or  not  the  electric  fields  pro- 
duced will   exceed  the  dielectric  breakdown  strength.     The  results  of 
CALCl   are  printed  by  subroutine  WRITER. 

For  the  graphical   output  of  the  results,  subroutine  OUTPUT  coordina- 
tes the  various  plotting  subroutines  required.     The  remainder  of  the 
required  data  cards  are  read  in  by  subroutines  OUTPUT  under  NAMELIST 
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control.  The  first  data  card  is  &TABL,  starting  in  column  2.  The 
variables  in  NAMELIST  TABL  are  read  into  two  vectors,  TABLE  and  LTG, 
in  the  format: 

"Array  name"  =  a  set  of  FORTRAN  constants,  separated  by  a  comma. 

The  format  for  the  elements  of  TABLE  is  TABLE(  )  =  ' ',  where  (-) 

indicates  a  card  column  available  for  a  data  item.  This  format  prepares 
the  data  to  be  used  in  labeling  the  graphical  output.  The  variables 
read  in  are: 

NAME  TYPE      DESCRIPTION 

Nozzle  radius  in  centimeters  ^ 
Initial  charge  density  in  C/m 
Maximum  fixed  potential  (GAPPOT) 
Dielectric  premittivity 
Breakdown  field  strength 
of  dielectric 
R*8      Characteristic  length  modeled 


TABLE  (1) 

R*8 

TABLE  (2) 

R*8 

TABLE  (3) 

R*8 

TABLE  (4) 

R*8 

TABLE  (5) 

.  R*8 

TABLE  (6) 


The  elements  of  the  vector  LTG  can  be  used  to  construct  a  grid  for  the 
plotted  output,  label  it  and  the  plotted  contours.  The  vector  LTG  is 
typed  LOGICAL  *  1.  The  variables  read  into  the  vector  LTG  are: 

LTG  (1)  1*1       .TRUE.  All  exterior  contour  segments 

and  those  interior  contour 
setments  which  represent  a 
local  maximum  will  be  labeled, 
•FALSE.  Omits  labeling  option. 

LTG  (2)  L*l       .TRUE.  Request  tic  marks  and  corres- 

ponding scale  values  one  inch 
apart  on  the  exterior  frame 
of  the  contour  graph. 
.FALSE.  Omits  "tic"  option. 

LTG  (3)  L*l       .TRUE.  Request  a  one  inch  by  one 

inch  straight  line  grid  to 
superimposed  on  the  contour 
graph . 
.FALSE.  Omits  the  grid  option. 

Again,  the  last' card  signals  the  end  of  the  NAMELIST.  The  format  is: 

&END 

with  the  symbol  ampersand  in  column  2  of  the  data  card. 
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APPENDIX  D:  MOBILITY  MEASUREMENTS 

The  apparatus  used  to  investigate  the  mobility  of  the  charged  particles 
is  shown  in  Figure  (19).  This  consists  of  a  fine  Monel  steel  grid  con- 
centric with  a  stainless  steel  cylinder.  The  EGD  process  under  investiga- 
tion involves  the  transport  of  negatively  charged  particles.  Hence  if 
a  voltage  is  applied  to  the  grid  and  varies  from  zero  to  some  value  below 
ground,  the  grid  will  repel  or  pass  negatively  charged  particles.  If  at 
the  same  time,  the  outer  cylinder  is  maintained  at  ground,  those  particles 
passed  by  the  grid  will  be  accelerated  across  the  gap,  giving  a  noticeable 
variation  in  current  flow  due  to  the  charge  arriving  at  the  outer  cylinder. 
From  the  definition  of  drift  velocity: 

V  =  k  E 

and  assuming  that  the  1  field  is  uniform: 

7  =  k  I  , 
where  1  is  the  spacing  between  the  grid  and  outer  cylinder.  But  since: 

^  -\ 

The  mobility  may  be  determined  by: 

^         tv 
By  varying  the  frequency  with  which  the  grid  voltage  changes  and  ensuring 
that  the  waveform  is  a  square  wave,  the  entire  range  of  mobility  possessed 
by  the  charged  particles  can  be  estimated,  assuming  that  mobility  is 
proportional  to  the  electric  field  between  the  cylinders.  This  is 
because  at  the  lower  frequencies,  all  particles  will  be  affected,  while 
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at  the  higher  frequencies,  or  small  t,  inertia  effects  will  deter  all 
but  the  smallest  or  most  highly  charged  particles  from  completing  their 
journey.  This  is  similar  to  the  method  of  Rutherford  as  outlined  by 
Darrow  [14]. 

The  measurements  leading  to  an  approximation  of  the  mobility  distri- 
bution were  made  at  a  reservoir  pressure  of  eight  psig,  over  a  wide  range 
of  temperatures,  238  -  296°F.  After  steam  flow  conditions  had  stabilzed, 
i.e.,  proper  pressure  and  constant  temperature,  the  device  was  inserted 
into  the  EGD  channel  to  such  a  position  that  it  was  axial  ly  symmetric 
with  the  steam  jet,  but  not  touching.  This  was  done  to  keep  condensation 
on  the  grid  to  a  minimum.  The  grid  was  powered  by  a  square  wave  genera- 
tor while  the  outer  ring  was  maintained  at  ground.  See  Figure  (18)  for 
a  schematic  of  the  setup.  Current  from  the  ring  was  fed  to  an  oscillo- 
scope for  observation.  It  was  found  that  the  ring  current  varied  con- 
siderably as  the  charged  particles  arrived  from  the  Monel  grid.  This 
variation  existed  over  a  range  of  fequencies.  At  high  temperatures, 
296°F,  the  range  of  current  excitation  appeared  to  be  concentrated  at 
a  single  point,  2.1  kHz.  As  the  temperature  of  the  steam  was  decreased, 
the  range  increased  considerably,  until  at  238°F,  the  current  excitation 
existed  from  an  upper  frequency  of  5  kHz  to  a  lower  frequency  of  .5  kHz. 

At  296°F,  the  mobility  corresponding  to  2.1  kHz  is  4.15  x  10'^ 

2  -3  2 

m  /V-s,  while  at  238° F,  the  mobility  range  runs  from  1  x  10   m  /V-s  to 

-4  2 
1  X  10   m  /V-s.  The  range  of  mobilities  appears  to  arise  from  both  the 

ability  of  large  particles  to  move  at  low  frequencies  and  the  lower  tem- 
perature producing  a  greater  distribution  in  particle  size.  However, 
further  work  remains  to  be  done  to  determine  both  size  distribution  and 
further  definition  of  the  mobility  distribution. 
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Figure  3.  Computer  Prediction  of  Parallel  Electrode  Normalized 
Potential  Distribution  at  Breakdown. 
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Figure  4.  Aerosol  Injector  and  Collector  Probe, 
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Pq  =  8  psig 


Pq  =  14  psig 


Figure  5.  Steam  Jet  Schlieren  Photographs 
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COMPUTER  PROGRAM 


C0MM0N/BL1/A(81,61 ) 

C0MM0N/BL2/B{81,61) 

C0MM0N/BL3/G(81t61) 

COMMON/BL^/NGLINEtNI, IW, IH 

C0MM0N/BL5/R(81) 

COMMON/BL6/NATRA0,RH0ZR0,PRMTVT,GAPPOTtBRKDWN,C0NVER, 
lNDZP0S,NGL-CS5,NR0WStNC0LST  DELU  t  NBLOCK,  NL  ,  NL  1 ,  HL 

COMMOfM/BLlO/NDI  MEM  ,  NDL  PNT  ,  NDLR  AD 

COMMON/BL11/XL1N(4000) ,YLIN(4000) 

NAMELIST/PARAM/NATRAr,RHCZRO,PRMTVT,GAPPOT, BRKDWN, 
lCONV£RtNOZPOb»NGLOSStNROWS,NCnLS,NDIMEN,DELU,NBLOCK, 
2NLtNLl ,NGL INE , NDLPNT , NDLRAD t I W , IH 

READ(5»PARAM) 

$ALPH^GAPP0T=)=100./(BRKDWN*C0NVER) 

$C0NST=(C0NVER/100.  )^=<^2*RH0ZR0/ (  PRMTVT*GAPPOT  ) 

WRITE (6tPARAM) 

HL  =  l./(  (NRCWS-l  )=^1.) 

WRITE(6,210)  $ALPHt«CONST 

DO  10  J=1,NC0LS 

DO  10  I=1,NRGWS 

A(I , J)=0.5 

B(I,J)=0.0 
10  G( It J)=0.0 
19  CALL  GSPRD($C0NSTta34) 

IF(NDIMEN.EQ.O)  GO  TO  34 

WRITE(6,240) 

WRITE(6,245)  NROWS 

WRITE(6,230)  {(J,R(J))  ,J  =  1, NROWS) 
34  CALL  GErjFIL(€,70) 

CALL  CALC1($ALPH) 

CALL  CUTPUT(&70) 
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GO  TO  70 

210 

FORMAT ( ' 

•  »2X, 

5F20.10) 

230 

FORMAT ( • 

•,T15 

,5CR(  ',13 

,' )=• ,P8.5,3X) 

»/) 

240 

FORMAT ( • 1 

•  ) 

245 

FORMAT ( 1 

•  tT49 

,'THE  JET 

RADIUS 

AT 

STATIONS 

1 

THRU   ', 

LI3,/) 

70' 

STOP 
END 
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SUBROUTINE  GSPRD  (  $CONST  ,i^) 
C0MM0N/DL3/G(bl,61) 
C0MMCr:/BL5/R{8  1  ) 

C0MM0N/BL6/NATRAD»RH0ZR0,PRMTVT,GAPP0TtBRKDWN,C0N\^ER, 
lNOZPOS,NGLCSS,NROWS,MCOLS,DELU,NBLOCK,NLtNLl,HL 
C0MM0N/6L10/NDIMEN,NDLPNT,NDLRAD 
Q=HL=<'-2*$C0NST 
GLOSS=0. 
IF(NDIMEN)    20,10,20 

10  DO    11    I=NOZPOS,NROWS 
DO    11    J=1,NCCLS 
IF(NGLOSS.EQ.l)     GLGSS=0. 

11  G(I,J>  =  (  l.-GL0SS)=^='!=2=*^Q 
IF(NOZPOS-l)     20,19,20 

19  RETURN    1 

20  IF(NOZPOS.EQ.NDLPNT)  GO  TO  24 
DO  21  I=1.NDLPNT 

21  R(I)=(NDLRAD-1)*HL 
K2=NDLPNT+1 

A=(  (NDLRAD-1  )=^HL)=<'-2 

B=(  (NATRAD-D'^HL  )=^-'^2 

VARA=1 A'B)/( (NDLPNT-NOZPOS  J*l. ) 

DO  23  I=K2,N0ZP0S 

23  R(I }=( VaRA*( I-N0LPNT)+A)**.5 
IF(NDIM5N.EQ.O)  RETURN  1 

24  A=(NR0WS-NGZP0S)=^-1. 
DO  35  I=NOZPOS,NROWS 
B=(I-NOZPOS)=i^l. 
X=B/A 

35  R(n  =  ,0445  +  ,0043-X  +  .C533*X**2 
B=R{NOZPOS) 
DO  50  I=NOZPOS,NROWS 
DO  40  J=1,NCDLS 
K1=J-1 
IF(R( i J-K1*HL)  43,42,42 

42  G(I,  J)  =  (B/R(I)-GL0SS)''^*2*Q 
GO  TO  40 

43  IF{R( 1 )-HL*(K1^.5) )  50,44,44 

44  G  (I  ,J)  =  (B/R(I  )-GLG,SS  )^=*2*Q 
GO  TO  50 

40  CONTINUE 

50  CONTINUE 
IF(NOZPOS.EQ.NDLPNT)  RETURN 
K1=NATRAD-1 
K3=2*NATRAD 

DO  90  I=K2,N0ZPCS 
DO  80  J=NDIRAD,K3 
IFCNGLOSS.EQ.l)  GLOSS=0. 
NRAD=J-NDLRAD 
IF(Kl-J)  60,51,51 

51  IF(R(  U-NRAD^^HL)  53,52,52 

52  Gd,  J)-(B/R(  I)-GL0SS)**2*Q 
GO  TO  60 

53  IF(R(I }-HL*(NRAD-.5) )  80,54,54 

54  G(  I,  J)  =  (B/R(  I)-GLOSS)**2'5=Q 
GO  TO  60 

60  NRAD-=2-NDLRAD-J 

IF{NRAD)  80,80,65 
65  G( I ,NRAO)=(B/R(I )-GL0SS)*^2*Q 
80  CONTINUE 
90  CONTINUE 

RETURN 

END 
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SUBROUTINE    GENFIL(=^) 

C0MM0M/BL1/A(8  1,61) 
C0.MM0N/BL2/B(81,61  ) 

C0iMM0N/BL3/G(  01,61) 
C0MM0N/3L4/N0LINE,Nlt IW, IH 

C0MM0N/BL6/NATRAD,RH0ZR0tPR'MTVTtGAPPGT,BRK0WN,C0NVER, 
1N0ZP0S,NGL0SS,MR0WS7N'C0LS,DELUtNBL0CK,NL,NL1tHL 
CG.MM0N/BL11/XLIM(4000)  ,YLIN(4000) 
XVAL  =  {  I  W^1.)/(NC0LS^^1.-1.) 
YVAL={  lH=:a.)/(NROWS-l.-l.) 
N1  =  0 


10 
20 

21 
70 
22 

200 


DO  23  I=1,NBL0CK 

READ(5,200)   Nil 

DO  10  Il=NIi,NI2 

DO  10  J1=NJ1,NJ2 

A{  I1,J1)  =  PV 

B(  II  ,J1)=PV 

G(  I1,J1)-GV 

IF(GV-98.98)  10,5,5 

N1=N1  +  1 

XLIN{K'l)=XyAL-(  Jl-1) 

YLIN(N1)=YVAL*( II- 1 ) 

CONTINUE 

CONTINUE 

IFtNOLINE)  22,21,22 

CALL  GUTPUT(&70) 

RETURN  1 

RETURN 

F0RMAT(4I3,2F12.3) 

END 


NI2,NJ1,NJ2,PV,GV 
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SUBROUTINE  CALC1($ALPH) 

C0MM0N/BL1/A(81,61) 
C0MM0N/bL2/B(8lT61) 

C0MM0N/BL3/G( 31,61) 
C0MM0N/BL5/R(81) 

COMMON /BL6/N AT  RAD tRH02R0,PRMTVT, GAP POT, BRK  OWN  tCONVER, 
lNOZPGS,hGLOSSTMRGWS,NCGLS,DELU,NBLOCK,NL,MLl,hL 
C0MM0N/BL10/NDIMEN,NDLPNT,NDLRAD 
ITCNT=1 
ITERS=1 
5  DO  20  J=1,NC0LS 
J1=J+1 
J2=J-1 
R1=J2-HL 
IF(Rl)  7,7,8 

8  ARG1=.5-HL/R1 

7  DO  15  I=1,NR0WS 
11=1+1 
12=1-1 

IF(G{ I,J).GT.98..AND.G(I,J )  .LT.IOO.)  GOTO  15 
IF( J.EQ.l.OR.J.EQ.NCCLS)  GO  TO  11 
IF( I .EQ.l.OR.I .EQ.NRCWS)  GO  TO  16 
IF(NDIMEM)  10,10,9 

9  A(I,  J)=B(I  ,J)  +  1.5*(.25"MG(  I  ,J)+ARG1*(A(  I,J1)-A(I,  J2)  ) 
1+A{I,J1)+A(I,J2)+A(I1,J)+A(I2,J))-B(I,J)) 

GO  TO  15 

10  A(  I,J)=B(I  ,J)  +  1.5*<.25*(G(  I,J)+A(I  ,  J1)+A(I,J2)+A{U,  J) 
1+A(I2, J) )-Bi  I, J)  ) 

GO  TO  15 

11  IF( J.EQ.NCOLS)  GO  TO  12 

A(I  ,  J)=B(I,l)+1.5^(1.333^A(I,2)-.3  33=f=A(  I,3)-B(I,1)  ) 
GO  TO  15 
IP  I  '^=  I  —  7 

A(I,  J)  =  B(I  ,J)+1.5*(1.333'i<A(  I ,  J  2  )-.  333'^  A(  I,J3)-B(I,J)) 

GO  TO  15 
16  IF( I .EQ.NROWS)  GO  TO  18 

A(1,J)=B(1,J)+1.5*(1.333^A{2,J )-.3  33*A(3, J )-B(l, J)  ) 

GO  TO  15 
18  13=1-2 

Ad,  J)  =  B(I,J)+1.5*(1.333*A(  I2,J)-.3  33-A(!3,J)-B(I,J)) 
15  CONTINUE 
20  CONTINUE 

DELMAX=0. 

DO  24  J=1,NC0LS 

DO  24  I=1,NR0WS 

DEL=A(  I,J)-B(I, J) 

ADEL=ABS(DEL) 

IF(ADEL.GT.DELMAX)  DELMAX=ADEL 
24  CONTINUE 

IF(DELU-DELMAX)  26,40,40 
26  DO  34  J=1,NCQLS 

DO  34  I=1,NR0WS 

34  B{ I, J)=A(I  ,Jj 

IF( ITCNT-100)  36,35,35 

35  WRITE(6,210)  ITERS, DELMAX 
ITCNT=1 
ITERS=ITERS+1 

GO  TO  5 

36  ITERS=ITERS+1 
ITCNT=ITCNT+1 
GO  TO  5   . 

40  WRITE(6,320)  ITERS 
WRITE(6,230) 
WRITE(6,232) 

CALL    WRITER( A,NROWS,NCGLS) 
ARG2  =  .  5/HL^"^$ALPH 
DO    45    J=1,NC0LS 
J1=J+1 
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41 
42 


43 
44 

45 


55 


109 

110 
113 

120 


49 
210 

230 
232 
234 
235 
237 
239 
320 

58 


J2=J-1 

J3=J+2 

J4=J-2 

DO  43  1=1 ,NROWS 

11=1+1 

12=1-1 

13=1+2 

14=1-2 

IF(I .EQ. 


1  = 
I, 


AZ=A( 12, 

GO  TO  42 
IFd.EQ. 
IF(I .EQ. 
IF(J.EQ. 
AR=A{ I, J 
GO  TO  44 
AR=0.0 
B(  I, J)=A 
G(  I,  J)=A 
CONTINUE 
WRITE(6, 
WRITE(6, 
CALL  WRI 
WRITE(6, 
WRITE(6, 
CALL  WRI 
EMAX=0.0 
DO  55  J= 
DO  5  5 
SUM=B( 
B(  I, J}=S 
IF(B(  I  ,J 
CONTINUE 
WRITE(6, 
WRITE(6t 
WRITE(6, 
CALL  WRI 
IF(NDIME 
DO  49  1= 
IF(R{I ) ) 
DO  113  J 
G(  I  , J)=l 
GO  TO  49 
XSLP  =  1  ./ 
DO  49  J= 
G( I, J)=X 
CONTINUE 
FORMAT ( • 
1F11.5,/) 
FORMAT ( 
FORiMAT( 
FORMAT ( 
FORMAT { 
FOR^'AT( 
FORMAT  ( 
FORMAT ( 


l.OR.I.EQ.NROWS) 
J)-A( II, J) 


GO  TO  41 


1)  AZ=A(3, J)+3.*A( 1,J)-4.*A(2, J) 
NROWS)  AZ  =  4.-A(  12 , J )-A ( 14, J )-3 . -A ( I , J  ) 

l.CR.J.EQ.NCOLS)  GO  TO  43 
2)-A( I,J1) 


RG2*AZ 
RG2-AR 


230) 

237) 

TER(B, NROWS, NCOLS) 

230) 

239) 

TERCG, NROWS, NCOLS) 

1, NCOLS 

1,  NROWS 

J)*^2  +  G(  I,  J)=^*2 

QRT(SUM) 

).GT.EMAX)  EMAX=B(I,J) 


230) 
234) 

235) 
TER( 
N)  5 
1,NR 
110 
=  1  ,  N 
.  +  1. 


EMAX 
B, NROWS, NCOLS) 
8, 5  8, 109 
OWS 

,  110,120 
COLS 
*( J-1) 


R{  I) 

1, NCOLS 

SLP=!=HL=<M  J-1) 

'  ,T5, 'AFTER'  , I5,1X,' ITERATIONS  —  DELMAX=» 


1'  ) 


,T49,'THE  MATRIX  OF  POTENTIALS',/) 
,T51,'THE  lEl  FIELD  MATRIX',/) 
,T42,' MAXIMUM  FIELD  STRENGTH  IS',F12,4,/) 
,T46, 'AXIAL  COMPONENTS  OF  |E|  FIELD',/) 
,T^7, 'RADIAL  COMPONENTS  OF  |E|  FIELD',/) 
,'THE  POTENTIAL  ITERATION  SCHEME  HAS', 


1T37,  'CONVERGED  IN', 15,'  ITERATIONS',/) 
RETURN 
END 
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SUBROUTINE  WRITER ( Z ,NROWS t NCOLS) 

DIMENSION  Z(8I,tl) 

Jl  =  l 

J2=10 
10  DO  15  I=1,NR0WS 

IF( J2,GE.NC0LS)  J2=NC0LS 
15  WRITE(6t220)  I , ( Z ( I , J ) , J=J 1 t J2 ) 

IF{ J2.EQ.NC0LS)  RETURN 

J1=J1+10 

J2=J2+10 

WRITE(6,230) 

GO  TO  10 
220  FORMATC  ',  13,  10^12.4) 
230  FORMATC  1*  ) 

END 


SUBROUTINE  CLVL ( CM , CLM ,NUML) 

C0MM0N/BL6/NATRADrRHCZR0,PRMTVT,GAPP0T,BRKDWN,C0NVER, 
lNOZPOS,NbLOSS,NRCWS,NCOLS,PELU,NBLOCK,NL,NLl,HL 

DIMENSION  CM(81,61 ),CLM(NUML} 

CMIN=CM(ltl) 

CMAX=CMIN 

DO  5  J=1tNC0LS 

DO  5  I=ltNROWS 

IF(CM(  I,  J)  .LT.CMIN)  CMIN=CM(I , J) 

IF(CM(  I ,J) .GT.CMAX)  CMAX  =  CM( I, J  ) 
i  CONTINUE 

NOW  DETERMINE  CONTOUR  LEVELS  TO  BE  PLOTTED. 

J=NUML 

I=NUML-1 

ANL=I*1. 

PLINT={CMAX-CMIN)/ANL 

NOW    FILL    THE    CONTOUR    LEVEL    VECTOR 

DC    6     I=ltJ 
.    CLMd  )  =CMIN+(I-1)*PLINT 

RETURN 

END 


SUBROUTINE  F LOP ( Z, NROWS » NCOLS ) 

DIMENSION  Z(81,61) 

IIVRT=NR0WS/2 

DO    3     I=1,IIVRT 

M=NRCWS-(I-1) 

DO    3    J=1tNC0LS 

SAVE=Z(I t J) 

Z( I, J)=Z(M,J) 

Z(M, J)=SAVE 

RETURN 

END 
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SUBROUTINE  OUT 
REAL=!=8  TITLEO 

1« 

2» 

3» 

4»GAPP0T 

5* 

6' 

7'AL  MAP 
REAL=!^8  TITLEA( 
REAL  *8  TITLEB 
REAL*8  TABLE(6 

.  LOGICAL-1  LTG( 
COMMGN/BLl/A 
C0MMGN/BL2/B(8 
C0MM0N/BL3/G 
C0MM0N/BL4/N0L 
C0MM0N/BL5/R{8 
C0MM0rvJ/BL6/NAT 

INOZPOS ,NGLGSS, 
C0MMGN/BL7/CL{ 
CGM1MON/BL8/CLI 
C0MMCN/BL9/CL2 
NAMELI ST/TABL/ 
READ(5,TABL) 
FIRST  FIGURE  C 
CALL  CLVL(AtCL 
NOW  FIGURE  CON 
CALL  CLVL(B,CL 
THE  CONTOUR  LE 
CL2( 1)=1.0 
WRITE(6,210) 
DO  3  I=l7NL 
3  WRITE(6,215)  I 
CALL  FLOP  (A,N 
CALL  FLOP  (3,N 
CALL  FLOP  (G,N 
TITLE(8) =TABLE 
TITLE! 11)=TABL 
TITLE( 14)=TASL 
TITLE( 17)=TA3L 
TITLE! 20)=TABL 
TITLE(23)=TABL 
IF(NOLINE)  7,6 

6  TITLE( 27)=TITL 
TITLE(28)=TITL 

7  CALL    CONTUR(A, 
11»  &20) 

NOW    PLOT     iEl     F 
TITLE( 27)=TITL 
TITLE(28)=TITL 
CALL    CONTURCB, 
lCL2t  It  £120) 
210    FORMAT (•     •  ,T26 
I'THE     lEl     FIELD 
215    FORMAT! '     • ,T33 
1F8.4,/) 
20    IF(NOLINE)     12t 
10    RETURN     1 
12    RETURN 
END 


PUT(*) 
0>/' 


RHOZRO 


CM.  BOH 


•  T» 


•  ,« 


El 


•PRMTVT 


•LENGTH 

•LEY 

•   BOX  89 

FIEL'  ,^D 


•NOZZLE 


•BRKDWN 


•  POTENTI 
/ 
MAP    •/ 
GEOMET^ , •RY  PLOT  •/ 


2)  /' 
(2)/ 
J 

3) 

(81, 

1,61 

(81, 

INE, 

1) 

RAD, 

NROW 

30) 

(30) 

(1  ) 

TABLE, LTG 


CNTOUR  LEVELS  TO 

,NL) 

TOUR  LEVELS  FOR  | 

1,NL1) 

VEL  VALUE  FOR  THE 


61) 

) 

61) 

Nl,  IW,  IH 

RHCZRO, PRMTVT 
S,NCGLS,DELU, 


,GAPPOT,BRKDWN,CONVER, 
NBL0CK,NL,ML1,HL 


BE  PLOTTED  FOR  MATRIX  A 
El  FIELD. 
CLOUD  BOUNDARY  IS  1.0 


,CL(  I ),I  ,CL1(I) 

ROWS, N COLS) 

ROWS,NCGLS) 

RGWS,NCOLS) 

(1) 

F(4) 

E(2) 

E(3) 

F(  5) 

E(6) 

,7 

EB(  1) 

EB{2) 

NRC  W S ,  NC  OL S  ,  8  1 ,  C  L  ,  NL  ,  T I T  L  E ,  LT G ,  G ,  CL  2  , 

I  ELD  CONTOURS. 

EA(1) 

EA(2) 

NR0VJS,NC0LS,81,CL1,NL1,TITLE,LTG,G, 


•THE  POTENTIAL  LEVELS  PLOTTED •,  T80 , 
PLOTTED^ ,//) 

)=^,F8.4,T87,'CL1(^,I2,^)=', 


LEVELS 
'CL( 


•,I2,^ 


10,12 
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SUBROUTINE  CONTUR  (  AM,  M  ,N,MX  ,CL  ,  NL,  T  ITL  E,  LTG,  G»  CL2  ,iML2, 

REAL'^e  TITLE(l) 

REAL^a  WIDTH/' WIDTH' /t HEIGHT/ 'HEIGHT'/, WHICH 

DIMENSION  AM(MX,1)  ,CL( 1) 

DIMENSION  GIMX,1),CL2( 1) 

DIMENSION  REC(900),   X(1800),  Y(1800) 

DIMENSION  IPT(3,3) ,INX(8) , INY( 8) 

COMMON  /OAYHOF/  MT , NT , N I , I  X, I Y , I DX , I DY , I SS , I T , I  V , NP , 
INQ, JT,PY,REC,CV, IPT, INX, INY,DL,RA,THE 

COMMON  /INTFAC/  X,Y 

DIMENSION  DITSX(5) ,DITSY{5 ) 

LOGICAL-1  LTG{ 1) , MINUS, LABL 

COMMON/TABL/  TABC ( 20 , 6 ) , JC 

COMMON/DITS/XMIN,YMIN,SLOPEX,SLOPEY, DI TSDX, DITSDY , 
1IDIR,LABL, MINUS 

C0MM0N/BL4/N0LINE,N1,IW, IH 

C0MM0N/BL10/NDIMEN,NDLPNT,NDLRAD 

C0MM0N/BL13/NSCN 

JC  =  0 

LABL=LTG(1) 

CHECK  IW  PARAMETER 

WHICH=WIDTH 


IFCIW)  1,1,2 

1 

WRITE! 6,60)  WHICH 

60 

FORMAT ( '0' ,T7, A8, ' 

OF  CONTOUR 

GRAPH 

ILLEGAL 

71 

WRITE (6,6^) 

64 

FORMAT! '0' ,T7, 'NO 

GRAPH 

WILL 

BE 

PRODUCED. 

1 

RETURN 

CHECK  IF  IW  IS  TOC 

)  WIDE 

2 

IF{IW-9)  3,3,40 

40 

WRITE (6,61) 

61 

FORMAT ( ' 0' ,T7, ' I W 

PARAMETER  i 

GREATER 

THAN 

9 

] 

L  WILL  SET  IW=9.' ) 

IW=9 

NOW  CHECK  IH  PARAMETER 

3 

IF(IH)  4,4,5 

4 

WHICH=HEIGHT 

CONTUR 


GO  TO  1 
5  01TSDX=(N-1.0)/IW 
0ITSDY=(-1.0+M)/IH 
XMIN=1.0 
YMIN=-M 

SL0PEX=1.0/DITSDX 
SL0PEY=1.0/DITSDY 
DITSX{ 1)=1.0 
DITSX(4)=1.0 
DITSX( 5)=1.0 
DITSX(2)=N 
DITSX(3)=N 
DITSY( 1)=-1.0 
DITSY(2)=-1.0 
DITSY( 5)=-1.0 
DITSY( 3)=-M 
DITSY(4)=-M 
DO  2011  1=1,5 

DITSX{  I  )  =  SLOPEX-"(DITSX(I  )-XMIN) 
2011  DITSY(  I  )=SLOPEY=^(DITSY(I)-YMIN) 
STARTP=(9.0-IW)/2.0 
CALL  PLOTS 

CALL  PL0T(STARTP,0.0,-3) 
CALL  LINE(DITSX,DITSY,5,1, 1)  . 
DITSX{  l)=DITSX(l)-.5 
DITSX(5)=DITSX(1 ) 
DITSX(4)=DITSX(4)-.5 
DITSX(2)=DITSX(2)+.5 
DITSX(3)=DITSX{3)+.5 
DITSY(  1)=DITSY{ l)  +  .5 


72 


DITSY(5)=0ITSY(1) 

DITSY(2)=DITSY(2)+.5 

DITSY( 3)=DITSY(3  )-.5 

DITSY(4)=DITSY(4)-.5 

CALL  LIME(DITSXtDITSY,5,1,1) 

SL0PEX=1.0/DITSDX 

SL0PEY=1.0/DITS0Y 

IENDX=SL0PEX^-=N+1 

IENDY=SL0PEY*M+1 

IF(.N0T.LTG{2) )  GO  TO  34 
C      DRAW  TIC  MARKS  ON  OUTER  FRAME 
C      START  ON  LEFT  EDGE  GOING  DOWNWARD 

IFLAG=0 

ZINGX=-.l 

ZINGY=0.0 

ZX=0.0 

ZY=-1.0 

CX=DITSX(1) 

CY=DITSY(l)-.5 

IEND=IENDY 
2222  IFLAG=IFLAG+1 

DO  2022  I=1,IEND 

CALL  PL0T(CX,CY,3) 

COORDX=CX+ZINGX 

COORDY=CY+ZINGY 

CALL  PL0T{C00RDX,C00RDY,2) 

CX=CX+ZX 
2022  CY=CY+ZY 

GO  TO  (21t22,23,24), IFLAG 
C      NOW  DO  THE  RIGHT  EDGE  GOING  DOWNWARD 

21  ZINGX=.l 
CX=DITSX{2) 
CY=DITSY(2)-.5 
GO  TO  2222 

C      NOW  DO  TOP  EDGE 

22  ZINGX=0.0 
ZINGY=.l 
ZX=1.0 
ZY=0.0 

CX=DITSX( l)+.5 
CY=DITSY(1) 
IEND=I ENDX 

GO  TO  2222 
C      NOW  DG  THE  BOTTOM  EDGE 

23  ZINGY=-.l 
CX=DITSX(4)+.5 
CY=DITSY(4) 
ZINGY=-.l 

GO  TO  2222 
C      NOW  LABEL  TIC  MARKS 
C      DO  X-DIRECTION  FIRST,  TOP  EDGE 
C      POSITION  PEN 

24  DELTAX=DITSDX 
IFLAG=0 
ZX=1.0 
ZY=0.0 

CX=DITSX( l)+.35 
CY=DITSY{I)+.12 

3033  IFLAG=IFLAG+1 

XZERO=1.0 

DO  3333  1=1, lEND 

CALL  NUMBER (CX,CY, . 14, XZERO ,0. 0 , 1 ) 

CX=CX+ZX 

CY=CY+ZY 
3333  XZERO=XZERO+DELTAX 

GO  TO  (31,32,33,34) , IFLAG 
C      LABEL  BOTTOM  EDGE  TIC  MARKS 

31  CX=DITSX(4)+.35 
CY=DITSY(4)-.19 
GO  TO  3033 

C      LABEL  LEFT  EDGE  OF  TIC  MARKS 

32  CX=DITSX(4)-.4 
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CY=DITSY(4)+.46 

DELTAX=DITSDY 
IEND=I ENDY 
ZX=0.0 
ZY=1.0 
GO  TO  3033 
;      NOW  LABEL  RIGHT  EDGE  TIC  MARKS 

33  CX=DITSX(3)+.12 
CY=DITSY(3)+.46 
GO  TO  3033 

;      CHECK  IF  GRID  DESIRED 

34  CALL  RESTOF(LTG, lENDXt IENDYTNL,AM,M,N,MX,CL,NL2,3,CL2t 
1STARTP,TITLEtDITSX,DITSY) 

IF(NOLINE)  2502,2505,2502 
2502  RETURN 
2505  RETURN  1 

END 


SUBROUTINE  REST0P(LTG, IENDX, IENDY, NL , AM » M,N ,MX ,CL , NL2 , 
2,TITLE,DITSX,DITSY) 

REAL'^^e  TITLE(l) 

DIMENSION  AM(MX,1) ,CL( 1) 

DIMENSION  G(MX,1) ,CL2(1) 

DIMENSION  DITSX{5) ,DITSY(5) 

LOGICAL^l    LTGd  )  tMINUS,LABL 

COMMON/TA6L/TABC(20»6)  ,JC 

C0Mi10N/DITS/XMI(\i,YMIN,SL0PEX,SL0PFYtDI  TSDXtDITSDY, 
lIDIR,LABL,niNUS 

C0MM0N/8L4/NaLINE,Nl,IW,Ui 

C0MM0M/BL10/NDIMEN,riDLPNT,NDLRAD 

COMMON/BL11/XLIN(4000) ,YLIN(4000) 

COMMON /BL13/NSCN 

IF( .N0T.LTG{3) )  GO  TO  35 
;      DRAW  INCH  BY  INCH  GRID 

IEND=IENDX-2 
;      POSITION  PEN 

IFLAG=0 

CX=DITSX( l)+.5 

CY=DITSY(l)-.5 

COORDX=0.0 

COORDY=-IH 

DX=1.0 

DY=0.0 
4044  DO  4444  I=ltIEND 

CX=CX+DX 

CY=CY+DY 

CALL  PL0T(CX,CYt3) 

ZX=CX+COORDX 

ZY=CY+COORDY 
4444    CALL    PLOT ( ZX , ZY , 2 ) 

IFdFLAG)    35,42,35 
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42 


35 
2500 

11 

20 


2002 
2004 
2003 


777 
778 


779 


IFLAG= 

IEND=I 

CY=DIT 

CX=DIT 

COORDX 

COORDY 

DX=0.0 

DY=1.0 

GO  TO 

CONTIN 

IF(NOL 

CALL  L 

IF(NOL 

DO  20 

CALL  S 

NSCN=1 

I F  (  N  D  I 

DO  200 

CALL  S 

IF( .NO 

IF( JC. 

DO  777 

COORDX 

COORDY 

CLEV=T 

CALL  N 

CONTIN 

CALL  S 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

RETURN 

CALL  S 

CALL  P 

CALL  P 

RETURN 

END 


ENDY-2 

SY(4)+.5 

5X(4)+.5 

=  IW 
=0.0 


4044 

UE 

INE)  11,2500,2500 

INE(XLIN,YLIN,N1 ,1 t-6) 

INE)     11,779,11 

I=1,NL 

CAN(AM,M,N,MX,CL(  I) ) 


MEN 
4  I 
CAN 
T.L 
EQ. 
1  = 
=  TA 
=  TA 
ABC 
UMB 
UE 
YME 
YM8 
YM6 
YMB 
YMB 
LOT 
LOT 


)    20( 

=  1,NI 

( G ,  M 

ABL) 

0)    G! 

1,JC 

BC(  I 

BCd 

(1,6 

ER{COORDX,CCORDY,.07,CLEV,0.0,3) 


102,2003,2002 
1L2 
,N,MX,CL2(I)) 

GO    TO    778 
;0    TO    778 

,4) 
»5) 
) 


OL(-STARTP, IH+1.0,.21,TITLF(25),0.0,48) 
0L( -START P, IH+1.5,.21,TITLE(19) ,3.0,48) 
0L(-STARTP,IH+2.0,.21,TITLE(13),0.0,48) 
OL(-STARTP, IH+2.5,  .21,TITLE(7),0.0,48) 
OL(-STARTP,IH+3.0 ,.21,TITLE(1),0.0,48) 
(-STARTP, IH+6.5,-3) 
E 


YMB0L(-STARTP,IH+1.0,.21,TITLE(2  5) ,0.0,48) 

LOT(-STARTP, IH+4.5,-3) 

LCTE 
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SUBROUTINE  SC AM { AM , M, N ,MX, CL J 


55 
57 


110 


15 
17 


DIMENSION  AM (MX 
DIMENSION  IPT(3 
COMMON  /DAYHGF/ 

COMMON  /INTFAC/ 
LOGICAL  -1  LABL 


1.570796 


t  ISS, ITtIV,NP,NQ 
INY,DL,RA,THE 


I     )tREC(900),       X{1800),    YdSOO) 

3) ,INX(8) , INY( 8) 

MT,NT,NI,  IX, lY, IDX,  IDY 

PY,REC,CV,  IPT,INX, 

X,Y 

MINUS 

ccmmon/dits/xmin,ymin,slopex,slopey,ditsdx,ditsdy, 
iidir, labl, minus 
c0mm0n/bl13/nscn 

D=0. 
R=l. 
TH    = 

NP  =  0 

DL  =  D 

RA  =  R 

THE=TH 

MT  =  N 

NT  =  M 

CV  =  CL 

IF(NSCN.EQ.l) . IZW=0 

IF(IZW-120631)     1,3,1 

IPT(1,1)=8 

2)=1 

3)=2 

1)=7 

3)=3 

1)=6 

2)=5 

3)  =4 

=-1 

=-1 

=  0 

=  1 

=  1 

=  1 

=  0 

=-1 
INY( 1)=0 
I  NY  {  2  )  =  1 
INY(3) =+1 
INY(4)=+1 
INY(5)=0 
INY(6)=-1 
INY{7)=-1 
INY(8)=-1 
IZW=120631 
3    XT=MT 

DO    58    J=l,900 
58    REC(J)=0 
ISS  =  0 
2    MT1=MT-1 
IDIR=1 
DO    110 
I  F  (  A  M  (  1 
IF(AM(  1 


IPT( 

IPT( 

IPT( 

IPT( 

IPT( 

IPT( 

IPT( 

INXd 

INX(2 

INX(3 

I  NX  (4 

INX(5 

INX{6 

INX(7 

INX(8 


IX=I+1 
IY  =  1 

IDX=-1 
IDY  =  0 

CALL  TRACE 
CONTINUE 
NT1=NT-1 
IDIR=2 
DO  20 
IF(AM( 
IF (AMI 
IX  =  MT 
IY=I+1 
IDX  =  0 


1=1, MTl 

,I)-CV)     55,110,110 

,I+1)-CV)     110,57,57 


(AM, MX) 


1=1, NTl 

I,MT)-CV)     15,20,20 

I+1,MT)-CV)    20,17,17 
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20 
22 


25 
27 


30 


35 
37 


40 


5 

7 

12 

9 
11 


10 


IDY  = 
CALL 

CONT 
ID  I  P. 
DO  3 
MT2  = 
IF(A 
IF(A 
IX  =  M 
IY  =  N 
IDX  = 
IDY= 
CALL 
CONT 
IDIR 
DO  4 
NT2  = 
IF(A 
IF(A 
IX  =  1 
IY  =  N 
IDX  = 
IDY  = 
CALL 
CONT 
IDIR 
ISS= 
NT1  = 
MTi  = 
DO  1 
DO  1 
IF(A 
IF(A 
COM  = 
IF  { 
DO  9 
IF  ( 
CONT 
IX  = 
IY  =  J 
IDX= 
IDY  = 
CALL 
CONT 
RETU 
END 


-I 

T 
IN 
=  3 
0 
MT 

r;( 

M( 

T2 

T 

1 

0 

T 
IN 
=  4 
0 

NT 
M( 
M( 

T2 

0 

1 

T 
IN 
=  5 
1 

NT 
MT 
0 
0 
M( 
M( 
IC 
NP 

I 
RE 
IN 
1  + 


RACE 
UE 


(AM, MX) 


I=1tMT1 

+  1-1 

NT,MT2)-CV)  25,30,30 

NT,MT2-1 )-CV)  30,27,27 

-1 


RACE 
UE 


(AM, MX) 


I=1,NT1 

+  1-1 

NT2,1)-CV)  35,40,40 

NT2-L,1)-CV)  40,37,37 


-1 


RACE 
UE 


(AM, MX) 


-1 
-1 
J  = 
1  = 
Jt 

J, 

0- 

J 

D= 
C( 
UE 
1 


2,NT1 

1,MT1 

I  }-CV)5, 

I+1)-CV) 

(I+l)+J 

12,11,12 

IjNP 
ID)-COM) 


10,10 
10,7, 


9,10,9 


-1 
0 

TRACE 
INUE 
RN 


(AM, MX) 
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SUBROUTINE  TRACE  (AM, MY) 

DIMENSJON  AM(MY,1  ),PEC(900),   X(.1800)»  YdSOO) 
DIMENSION  IPT(3,3),  INX( 8) , INY ( 8  ) 

COMMON  /DAYHGF/  MT  ,NT  ,  NI  »  I  X  ,  I Y  ,  I  DX  ,  I  DY  ,  I  SS  ,  IT  ,  I\/ ,  NP  , 
IN  tJT,PYtREC,CV,IPT,INX,INY,DL,RA,THE 
COMMON  /INTFAC/  X,Y 
PY=0,0 

RC=   COS  (THE)=^-RA 
RS=   SIN  (THE)=;=RA 
501  JT=0 
N=0 
IXO=IX 
IYO=IY 
ISX=IDX+2 
ISY=IDY+2 
IS=IPT(ISX,ISY) 
JTB=0 
ISO=IS 
IF(IS0-8)18,18,17 

17  IS0=IS0-8 

18  IT=0 

5  CONTINUE 

CALL  CALC   (AM, MY) 

NZ  =  N 

N=NZ 

IF  (IT+JT-1)  49,49,47 
47  XS=X(N-1) 

YS=Y(N-1) 

X(N-1)=X(N) 

Y(N-1)=Y(N) 

X(N) =XS 

Y(N)=YS 
49  IS=IS+1 

JT  =  IT 
9  IF  (IS-9)  8,7,7 

7  IS=IS-8 

8  IDX=1NX(IS) 
IDY=INY( IS) 
IX2=IX+1DX 
IY2=IY+IDY 
JTB=JTB+1 

IF  (JTB-1799)  51,51,308 
308  PRINT  103,CV,X(N) ,Y(N) 

103  F0RMAT(1H0,23HA  CONTCUR  LIME  AT  LEVEL,E12.5, 
121H  WAS  TERMINATED  AT  X=,E12.5,3H  Y=,E12.5/ 
2  40H  BECAUSE  IT  CONTAINED  MORE  THAN  1799  PLOT  POINTS 
RETURN 
51  CONTINUE 

IF  (ISS)  10,10,20 

20  IF (TX- I XO)  12,21,12 

21  IF(IY-IYO)  12,22,12 

22  IF(IS-ISO)  12,23,12 

23  CONTINUE 

CALL  CALC   (AM, MY) 
GO  TO  73 

10  IF(IX2)  13,50,13 

13  IF  (IX2-MT)  19,19,50 

19  IF  ( IY2)  11,50,11 

11  IF  (IY2-NT)  12,12,50 

12  IF(CV-AM( IY2,IX2) )  206,206,5 
206  IF  {  IDX^=-2  +  IDY''=^!'2-l)  213,6,213 

213  DCP=(AM{ IY,IX)+AM( lY, IX2)+AM( IY2,IX)+AM( I Y2 , I X2) ) /4. 0 
IF  (DCP-CV)  5,217,217 

217  IF  (INXdS-D)  214,215,214 

214  IX=IX+IDX 
IDX=-IDX 
PY=2.0 

CALL  CALC   (AM, MY) 
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215 


6 
306 

16 


50 

307 

73 

74 


IX=IX+IDX 

GO    TO    6 

IY=I Y+IDY 

IDY=-IDY 

PY=2.0 

CALL    CALC       (AMtMY) 

IY=I Y+IDY 

IF(AM( lY, IX-1)-CV)     306,16tl6 

NP=NP+1 

REC(NP)=100*IX+IY 

IS=IS+5 

IX=I X2 

IY=IY2 

GO    TO    9 

XT=MT 

IF(AM{ lY, IX-1)-CV)     307,73,73 

NP=rjp+i 

REC(NP)=100*IX+IY 
DO    74     I=1,N 
X(I)=X(I  )  +  RC*Y(I) 
Y(  I)=RS=^Y{I) 
CALL    PLOTT(NtCV) 
RETURN 
END 
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SUBROUTINE    C  ALC  (  AiM  t  MY  ) 

DIMENSICM  A.M(MY,1  )tREC(900)»   X(1800)t  Y(1800) 
DIMENSION  IPT(3t3) , INX(8), INY(8) 

COMMON  /DAYHOF/  MT , NT , N I , I  X , I Y ,  I DX » I DY , ISS ,  IT , W , NP , N, 
lJT,PY,REC,CV,IPT,INXTlNY,DLtRA,THE 
COMMON  /INTFAC/  X,Y 
IT  =  0 
N=N+1 
IF  {IDX**2  +  IDY^^=!=2  -1)  20,1t20 

1  IF  (IDX)  10,2,10 

2  X(N)=IX 
Z=IY 

. IY2=I Y+IDY 

DY=IOY 
41  Y(N)  =  (  (AM{IY,IX)-CV)/(  AM(IY,IX)-AM(  IY2,  IX)  )  )=^DY  +  Z 

RETURN 
10  Y(N)=IY 

W=IX 

DX=IDX 

IX2=IX+I0X 
44  X(N)=( (AM( IY,IX)-CV)/( AM(IY,IX)-AM( lY, 1X2) ) )*DX+W 

RETURN 

20  IX2=IX+IDX 
IY2=IY+IDY 
W=IX 

Z=IY 

DX=IDX 

DY=IDY 

DCP=(  AM(  IY,IX)+AM{  lY,  IX2)+AM(  I  Y2  , 1  X  ) +AM(  I Y2  ,  I  X2  )  }  /4.  0 

IF    (PY-2.0)    24,21,24 

24  IF    (DCP-CV)    21,21,25 

21  AL  =  AM(  lY,  IX)-DCP 

23    V=.5*(AL+DCP    -CV)/AL 

27  X(N)=V-DX+W 
Y(N)=V-DY+Z 
PY=0.0 
RETURN 

25  IT=1 

AL=AM( IY2,IX2)-DCP 
33    V=.5=!^(  AL  +  DCP'-CV)/AL 

28  X{N)=-V'-DX+W  +  OX 
Y{N)=-V^DY+Z  +  DY 
Y(NJ=-V^=DY+Z  +  DY 
RETURN 

END 
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SUBROUTINE  PLOTT(NPtCV) 
COMMON/INTFAC/X( 1800) , Y( 1800) 
LOGICAL-1  .MINUS,  LABL 
COMMCN/TABL/  TABC ( 20  ,  6 )  ,  JC 

COMMON /D I TS/XM I N,YMIN, SLOP  EX, SLOPE Y,DITSDX,DITSDY, 
1IDIR,LABL, MINUS 
C  SCALE    POINTS    FOR    PLOT    ROUTINE 

DO    100     1  =  1, NP 
X(  I  )=SLCPEX^MX(I  )-XMIN) 
100    Y( I )=SLOPEY*(-Y( I)-YMIN) 
CALL    LINE(X,Y,NP,1, 1) 
IF{ .NOT.LABL)     RETURN 
DIR=0.0 
GO    TO     (1,2,3,4,6),     ICIR 

1  DIR=90. 

2  COCRDX=X( 1  ) 
C00RDY=Y(1  ) 

5  CALL  NUMBER(COORDX,CCaRDY, .07,CV,DIR,3) 
RETURN 

C      MOVE  PEN  DOWN  ONE  HALF  INCH 

3  DIR=90. 
C00RDX=X(1  ) 
C00RDY=Y(1  )-.3 
GO  TO  5 

C      MOVE  PEN  TO  THE  LEFT 

4  C00RDX=X(1  )-.3 
C00RDY=Y{1  > 

GO  TO  5 
C      SEARCH  FOR  XMAX , XM IN , YMAX, YMIN , AND  SAVE  YMINX 

6  XMAX=X(1) 
SMIN=XMAX 
YMINX=Y(1) 
YMAX=YMINX 
VMIN=YMINX 

DO  200  1=2, NP  • 
IF(X(I ) .GT.XMAX)  XMAX=X(I) 
IF(Y{  I  )  .LT.VMIN)  VMir!=Y(I) 
IF(Y( I ) .GT.YMAX)  YMAX=Y(I) 
1F{X( I ) .GE.SMIN)  GO  TO  200 
SMIN=X{I) 
YMINX=Y( I ) 
200  CONTINUE 
C      JC=NUMBER  OF  ENTRIES  IN  TABC 
IF(JC)  400,500,400 
400  DO  900  1  =  1, JC 

IF(XMAX.LT.TABC( 1,1)  .AND .YMAX . LT .TABC (  1,2)  .AND.VMIN. 
IGT.TABCd  ,3)  .AND.SMIN.GT.TABCd  ,4)  )  GO  TO  700 
900  CONTir:UE 
C      DID  NOT  FIND  THIS  CONTOUR  TO  BE  INTERIOR  TO  ANOTHER 
C      CHECK  IF  EXTERIOR 
DO  1000  1=1, JC 

IF(XMAX.GT.TABC( 1,1) .AND ./MAX. GT .TABC ( 1 , 2 )  . AND. VM IN. 
1LT.TABC(  I  ,3)  .AND.SMIN.LT.TABCd  ,4)  )  GO  TO  800 
1000  CONTINUE 
500  IF  (JC.EQ^20)  RETURN 
JC=JC+1 
MC  =  JC 
600  TABC  (MCI)  =XMAX 
TABC(MC,2)=YMAX 
TABC(MC,3)=VMIN 
TABC{MC,4) =SMIN 
TABC(MC, 5)=YMINX 
TABC(MC,6)=CV 
RETURN 
C      CHECK  IF  THIS  INTERIOR  ONE  IS  OF  HIGHER  LEVEL 
700  IF{CV.LE.TABC{ 1,6)  )  RETURN 
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2000  MC=I 

GO  TO  600 
;      CHECK  IF  LEVEL  OF  THIS  EXTERIOR 
800  IFCCV.LT.TABCd  ,6)  )  RETURN 

GO  TO  2000 

END 


ONE  IS  HIGHER 
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